home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATHEMAT / 1769.ZIP / PFSAI.FOR < prev    next >
Text File  |  1988-10-05  |  58KB  |  2,949 lines
  1. $include: 'picom.for'
  2.     open(unit=7,file='algin',status='old')
  3.     open(unit=8,file='algout',status='unknown')
  4.     call datim
  5.     itbeg=itime
  6.     do 4 i=1,nvar
  7. 4    isb(i)=0
  8.     ios(1)=1
  9.     iosi(1)=1
  10.     ipv(1)=1
  11.     pv(1,1)='?'
  12.     isb(1)=1
  13.     ise(1)=1
  14.     itb(1)=1
  15.     ite(1)=1
  16.     ivn(1)=2
  17.     ivp(1)=1
  18.     nd=0
  19.     nv=1
  20.     ns=1
  21.     np=2
  22.     nt=1
  23.     nb=7
  24.     maxe=9999
  25.     maxnt=0
  26.     maxite=0
  27.     itetop=1
  28.     ndump=0
  29.     minns=9999
  30.     lintot=0
  31.     idflg=1
  32.     call defh
  33.     ma=ntot
  34.     mb=nvar
  35.     mc=nterm
  36.     md=ninch
  37.     me=ninl
  38.     mg=nder
  39.     mh=nlab
  40.     mi=ndind
  41.     mj=mitay
  42.     write(8,1)ma,mb,mc,md,me,mg,mh,mi,mj
  43. 1    format(' version 85/6/20 ',
  44.      +   ' NTOT=',i6,' NVAR=',i4,' NTERM=',i5,' NINCH=',i5
  45.      +   ,/,' NINL=',i6,' NDER=',i3,' NLAB=',i3,' NDIND=',i2
  46.      +   ,' MITAY=',i2)
  47.     call monitr
  48.     call datim
  49.     write(8,90)nd,nv,ns,maxnt,maxite,lintot
  50. 90    format(' nu. derivatives=',i3,'  nu. variables=',i3,
  51.      +    '  nu. sums=',i3,/,' max nu. terms=',i5,
  52.      +  '  max total nu. variables=',i6,'  input lines=',i5)
  53.     t=(itime-itbeg)/100.
  54.     write(8,91)t
  55. 91    format(' time used=',f9.2,' seconds')
  56.     stop
  57.     end
  58.  
  59.  
  60.     subroutine addrat(n1,i1,n2,i2,n3,i3)
  61. c  n3/i3=n1/i1 + n2/i2
  62.     implicit integer*4 (i-n)
  63.     i=i1*i2
  64.     n=n1*i2+n2*i1
  65.     call gcdout(n,i)
  66.     n3=n
  67.     i3=i
  68.     return
  69.     end
  70.  
  71.  
  72.     subroutine adds(is,js,iin,iid,ians)
  73. $include: 'picom.for'
  74.     integer*4 in1,id1,in,id,inn,idd,insav,idsav
  75.     in=iin
  76.     id=iid
  77.     nt1=nt+1
  78.     ib=isb(is)
  79.     ie=ise(is)
  80.     jb=isb(js)
  81.     je=ise(js)
  82.     i=itb(ie)
  83.     if(ivn(i).ne.2 .or. ivp(i).le.maxe)goto 3
  84.     do 1 it=ib,ie
  85.     i=itb(it)
  86.     if(ivn(i).eq.2 .and. ivp(i).gt.maxe)goto 2
  87. 1    continue
  88. 2    ie=it-1
  89. 3    i=itb(je)
  90.     if(ivn(i).ne.2 .or.ivp(i).le.maxe)goto 7
  91.     do 4 it=jb,je
  92.     i=itb(it)
  93.     if(ivn(i).eq.2 .and. ivp(i).gt.maxe)goto 5
  94. 4    continue
  95. 5    je=it-1
  96. 7    if(ib.gt.ie)goto 70
  97.     if(jb.gt.je)goto 12
  98.     if(icomp(ib,jb))8,60,30
  99. 8    jl=jb
  100.     if(jb.eq.je)goto 10
  101.     if(icomp(ib,je))10,16,18
  102. 10    call copyt(jb,je,iin,iid)
  103. 12    call copyt(ib,ie,1,1)
  104. 14    call defsum(ians,nt1,nt)
  105.     return
  106. 16    insav=itn(je)
  107.     idsav=itd(je)
  108.     idd=id
  109.     if(in.lt.0)idd=-idd
  110.     inn=iabs(in)
  111.     call mulrat(itn(ib),itd(ib),idd,inn,in1,id1)
  112.     call addrat(in1,id1,itn(je),itd(je),itn(je),itd(je))
  113.     call copyt(jb,je,iin,iid)
  114.     itn(je)=insav
  115.     itd(je)=idsav
  116.     ib=ib+1
  117.     jb=je+1
  118.     goto 7
  119. 18    jr=je
  120. 20    if(jr-jl.eq.1)goto 28
  121.     jc=(jl+jr)/2
  122.     if(icomp(ib,jc))22,26,24
  123. 22    jl=jc
  124.     goto 20
  125. 24    jr=jc
  126.     goto 20
  127. 26    insav=itn(jc)
  128.     idsav=itd(jc)
  129.     idd=id
  130.     if(in.lt.0)idd=-idd
  131.     inn=iabs(in)
  132.     call mulrat(itn(ib),itd(ib),idd,inn,in1,id1)
  133.     call addrat(in1,id1,itn(jc),itd(jc),itn(jc),itd(jc))
  134.     call copyt(jb,jc,iin,iid)
  135.     itn(jc)=insav
  136.     itd(jc)=idsav
  137.     ib=ib+1
  138.     jb=jc+1
  139.     goto 7
  140. 28    call copyt(jb,jl,iin,iid)
  141.     jb=jr
  142. 30    il=ib
  143.     if(ib.eq.ie)goto 31
  144.     if(icomp(jb,ie))31,36,38
  145. 31    call copyt(ib,ie,1,1)
  146. 32    call copyt(jb,je,iin,iid)
  147.     goto 14
  148. 36    insav=itn(ie)
  149.     idsav=itd(ie)
  150.     call mulrat(in,id,itn(jb),itd(jb),in1,id1)
  151.     call addrat(in1,id1,itn(ie),itd(ie),itn(ie),itd(ie))
  152.     call copyt(ib,ie,1,1)
  153.     itn(ie)=insav
  154.     itd(ie)=idsav
  155.     jb=jb+1
  156.     ib=ie+1
  157.     goto 70
  158. 38    ir=ie
  159. 40    if(ir-il.eq.1)goto 48
  160.     ic=(il+ir)/2
  161.     if(icomp(jb,ic))42,46,44
  162. 42    il=ic
  163.     goto 40
  164. 44    ir=ic
  165.     goto 40
  166. 46    insav=itn(ic)
  167.     idsav=itd(ic)
  168.     call mulrat(in,id,itn(jb),itd(jb),in1,id1)
  169.     call addrat(in1,id1,itn(ic),itd(ic),itn(ic),itd(ic))
  170.     call copyt(ib,ic,1,1)
  171.     itn(ic)=insav
  172.     itd(ic)=idsav
  173.     jb=jb+1
  174.     ib=ic+1
  175.     goto 7
  176. 48    call copyt(ib,il,1,1)
  177.     ib=ir
  178.     goto 8
  179. 60    insav=itn(ib)
  180.     idsav=itd(ib)
  181.     call mulrat(in,id,itn(jb),itd(jb),in1,id1)
  182.     call addrat(in1,id1,itn(ib),itd(ib),itn(ib),itd(ib))
  183.     call copyt(ib,ib,1,1)
  184.     itn(ib)=insav
  185.     itd(ib)=idsav
  186.     ib=ib+1
  187.     jb=jb+1
  188.     goto 7
  189. 70    if(jb.gt.je)goto 14
  190.     goto 32
  191.     end
  192.  
  193.  
  194.     subroutine bomb(i,h)
  195. $include: 'picom.for'
  196.     character*8 h
  197.     goto(100,20,30,40,50),i
  198. 20    write(8,21)
  199. 21    format(' parameter error')
  200.     goto 100
  201. 30    write(8,31)h
  202. 31    format(1x,a8)
  203.     goto 100
  204. 40    write(8,41)h,h
  205. 41    format(1x,a8,' is too small for your problem.',
  206.      +       '  Recompile with bigger ',a8)
  207. 50    continue
  208. 100    write(8,101)line,ap
  209. 101    format(' BOMB: input line nu.',i5,' was:',/,1x,127a1)
  210.     write(8,90)nd,nv,ns,maxnt,maxite,line
  211. 90    format(' nu. derivatives=',i3,
  212.      +   '   nu. variables=',i3,' nu. sums=',i3,/,'  max nu. terms=',
  213.      +      i5,'  max total nu. variables=',i6,' input lines=',i5)
  214.     call datim
  215.     if(ndump.ne.0)call dump
  216.     t=(itime-itbeg)/100.
  217.     write(8,91)t
  218. 91    format(' time used=',f9.2,' seconds')
  219.     stop
  220.     end
  221.  
  222.  
  223.     subroutine bkup(j,k)
  224. $include: 'picom.for'
  225.     if(j.eq.k)goto 2
  226.     do 1 i=j,k-1
  227.     itb(i)=itb(i+1)
  228.     ite(i)=ite(i+1)
  229.     itn(i)=itn(i+1)
  230. 1    itd(i)=itd(i+1)
  231. 2    k=k-1
  232.     return
  233.     end
  234.  
  235.  
  236.     subroutine call
  237. $include: 'picom.for'
  238.     call nxt(3,nchpl,' ',j1)
  239.     call numl(j1,nchpl,j2,j3,j4)
  240.     if(linlab(j4).eq.-1)call bomb(3,'call    ')
  241.     linret=line
  242.     line=linlab(j4)
  243.     return
  244.     end
  245.  
  246.  
  247.     subroutine coef(is,jv,in,id)
  248. c is=start, on return jv points beyond coef
  249. $include: 'picom.for'
  250.     integer*4 in,id,ii,ij
  251.     id=1
  252.     isgn=1
  253.     call srchnb(is,nchpl,j1)
  254.     if(j1.ne.-1)goto 5
  255. 1    jv=-1
  256. 2    in=isgn
  257.     return
  258. 5    isgn=-1
  259.     if(a(j1).eq.'-')goto 30
  260.     isgn=1
  261.     if(a(j1).eq.'+')goto 30
  262.     j2=j1
  263. 7    jv=j2
  264.     call intval(j2,i)
  265.     if(i.eq.-1)goto 2
  266.     call nxtnnu(j2,nchpl,j3)
  267.     call getint(j2,j3-1,in2,in)
  268.     j5=j3
  269.     if(a(j3).eq.'.')goto 10
  270.     if(a(j3).eq.'/')goto 25
  271.     if(a(j3).eq.' ')goto 20
  272.     goto 15
  273. 10    call intval(j3+1,i)
  274.     if(i.ne.-1)goto 16
  275.     j5=j3+1
  276.     if(a(j5).eq.'/')goto 25
  277.     if(a(j5).eq.' ')goto 20
  278. 15    jv=j5
  279.     in=in*isgn
  280.     return
  281. 16    call nxtnnu(j3+1,nchpl,j7)
  282.     call getint(j3+1,j7-1,in2,ii)
  283.     ij=10**(j7-j3-1)
  284.     in=ij*in+ii
  285.     id=ij
  286.     call gcdout(in,id)
  287.     goto 26
  288. 20    call srchnb(j5,nchpl,j1)
  289.     j5=j1
  290.     if(j5.eq.-1)goto 15
  291.     if(a(j5).eq.'/')goto 25
  292.     goto 15
  293. 25    call nxtnb(j5+1,nchpl,j6)
  294.     call intval(j6,i)
  295.     if(i.eq.-1)goto 15
  296.     call nxtnnu(j6,nchpl,j7)
  297.     call getint(j6,j7-1,in2,id)
  298. 26    if(a(j7).eq.'.')j7=j7+1
  299.     call srchnb(j7,nchpl,jv)
  300.     in=in*isgn
  301.     return
  302. 30    call srchnb(j1+1,nchpl,j2)
  303.     if(j2.eq.-1)goto 1
  304.     goto 7
  305.     end
  306.  
  307.  
  308.     subroutine comput
  309. $include: 'picom.for'
  310. 300    call nxt(1,nchpl,' ',j1)
  311.     call numv(j1,nchpl,j2,j3,j4)
  312.     call nxt(j3,nchpl,'=',j5)
  313.     call srch(j5,nchpl,'/',j6)
  314.     if(j6.eq.-1)goto 320
  315. c  /
  316.     call nxtnb(j5+1,j6-1,j7)
  317.     if(a(1).eq.'d')goto 301
  318.     call numv(j7,j6-1,j8,j9,j10)
  319.     call numv(j6+1,nchpl,j11,j12,j13)
  320.     if(isb(j10).eq.0 .or. isb(j13).eq.0)call bomb(3,'not sum ')
  321.     call recip(j13,1)
  322.     nt1=nt+1
  323.     call mult(isb(j10),ise(j10),isb(1),ise(1))
  324.     call defsum(j4,nt1,nt)
  325.     call outqu(j4)
  326.     minns=min0(minns,ios(1))
  327.     it=isb(1)
  328.     ise(1)=it
  329.     ite(it)=itb(it)
  330.     k=itb(it)
  331.     ivn(k)=2
  332.     ivp(k)=0
  333.     call pack
  334.     return
  335. c  compute sj4=dsj10/dfj14
  336. 301    call numv(j7+1,j6-1,j8,j9,j10)
  337.     if(a(j7).ne.'d')call bomb(3,'comput1 ')
  338.     if(idflg.ne.0)call dertab
  339.     call nxt(j6,nchpl,'d',j11)
  340.     call numv(j11+1,nchpl,j12,j13,j14)
  341.     call srch(j13,nchpl,'=',j15)
  342.     iord=1
  343.     if(j15.eq.-1)goto 311
  344.     call nxtnb(j15+1,nchpl,j16)
  345.     call intval(j16,j17)
  346.     if(j17.ne.-1)goto 360
  347.     call numv(j16,nchpl,j19,j20,j18)
  348.     if(isb(j18).eq.0)call bomb(3,'comput14')
  349.     iord=itn(isb(j18))
  350.     if(itd(isb(j18)).ne.1)call bomb(3,'comput15')
  351.     if(iord.gt.0)goto 311
  352.     call bomb(3,'comput16')
  353. 360    call getexp(j15+1,iord,j16,0)
  354.     if(j16.ne.-1 .or. iord.le.0)call bomb(3,'comput2 ')
  355. 311    do 313 ii=1,iord
  356.     if(ii.eq.1)call deriv(j10,j14,1)
  357.     if(ii.ne.1)call deriv(1,j14,1)
  358.     call outqu(1)
  359.     if(itn(isb(1)).eq.0)goto 314
  360. 313    continue
  361. 314    call swap(j4)
  362.     return
  363. c  not /
  364. 320    call srch(j5,nchpl,'+',j6)
  365.     if(a(1).eq.'d')call bomb(3,'comput3 ')
  366.     id12=1
  367.     if(j6.ne.-1)goto 321
  368.     call srch(j5,nchpl,'-',j6)
  369.     in12=-1
  370.     if(j6.ne.-1)goto 322
  371.     call srch(j5,nchpl,'*',j6)
  372.     if(j6.ne.-1)goto 339
  373.     call srch(j5,nchpl,'^',j6)
  374.     if(j6.ne.-1)goto 339
  375. c  compute sj4=sj12
  376.     call numv(j5,nchpl,j6,j11,j12)
  377.     if(isb(j12).eq.0)call bomb(3,'comput4 ')
  378.     call copys(j12,j4)
  379.     call outqu(j4)
  380.     return
  381. c  compute sj4=sj9  + or -  sj12
  382. 321    in12=1
  383. 322    call numv(j5+1,j6-1,j7,j8,j9)
  384.     call numv(j6+1,nchpl,j10,j11,j12)
  385.     call adds(j9,j12,in12,id12,j4)
  386.     call outqu(j4)
  387.     return
  388. c  compute sj4=sj9*sj12
  389. 339    call numv(j5+1,j6-1,j7,j8,j9)
  390.     if(a(j6).eq.'^')goto 350
  391.     if(a(j6+1).eq.'*')goto 349
  392.     call numv(j6+1,nchpl,j10,j11,j12)
  393.     nt1=nt+1
  394.     call mult(isb(j9),ise(j9),isb(j12),ise(j12))
  395.     call defsum(j4,nt1,nt)
  396.     call outqu(j4)
  397.     return
  398. 349    j6=j6+1
  399. 350    call numv(j6+1,nchpl,j10,j11,j12)
  400.     call expnd(j9,j12,j4)
  401.     return
  402.     end
  403.  
  404.  
  405.     subroutine copyc(a,b,n)
  406.     integer*2 i,n
  407.     character*1 a,b
  408.     dimension a(2),b(2)
  409.     do 1 i=1,n
  410. 1    b(i)=a(i)
  411.     return
  412.     end
  413.  
  414.  
  415.     subroutine copys(j1,j2)
  416. $include: 'picom.for'
  417.     nt1=nt+1
  418.     call copyt(isb(j1),ise(j1),1,1)
  419.     call defsum(j2,nt1,nt)
  420.     return
  421.     end
  422.  
  423.  
  424.     subroutine count
  425. $include: 'picom.for'
  426.     call nxt(3,nchpl,' ',j1)
  427.     call numv(j1,nchpl,j3,j4,j5)
  428.     do 12 i=nchpl,j4+1,-1
  429.     if(a(i).ne.' ')goto 14
  430. 12    continue
  431.     call bomb(3,'count1  ')
  432. 14    do 16 j=i,j4+1,-1
  433.     if(a(j).eq.' ')goto 18
  434. 16    continue
  435. 18    call numv(j,i,j6,j7,j8)
  436.     if(isb(j8).eq.0)call bomb(3,'count2  ')
  437.     if(isb(j5).ne.0)goto 20
  438.     k=itetop+1
  439.     call incnt
  440.     itb(nt)=k
  441.     call nite(k)
  442.     itn(nt)=ise(j8)+1-isb(j8)
  443.     itd(nt)=1
  444.     ivn(k)=2
  445.     ivp(k)=0
  446.     call defsum(j5,nt,nt)
  447.     return
  448. 20    ia=isb(j5)
  449.     ib=ise(j5)
  450.     k=itb(ia)
  451.     l=ite(ia)
  452.     ivn(k)=2
  453.     ivp(k)=0
  454.     ise(j5)=ia
  455.     ite(ia)=k
  456.     itn(ia)=ise(j8)+1-isb(j8)
  457.     itd(ia)=1
  458.     if(ia.ne.ib .or. k.ne.l)minns=min0(minns,ios(j5))
  459.     call pack
  460.     return
  461.     end
  462.  
  463.  
  464.     subroutine copyt(ia,ib,in,id)
  465. c  copy terms ia-ib to nt+1,..., coefficients are multiplied by in/id
  466. $include: 'picom.for'
  467.     integer*4 in4,id4
  468.     in4=in
  469.     id4=id
  470.     if(ib.lt.ia)return
  471.     if(ia.lt.1 .or. ib.gt.nt)call bomb(3,'copyt2  ')
  472.     icop=0
  473.     if(id.eq.1 .and. in.eq.-1)icop=-1
  474.     if(id.eq.1 .and. in.eq. 1)icop= 1
  475.     do 15 i=ia,ib
  476.     ja=itb(i)
  477.     if(ivn(ja).eq.2 .and. ivp(ja).gt.maxe)goto 15
  478.     k=itetop
  479.     call incnt
  480.     itn(nt)=itn(i)*icop
  481.     itd(nt)=itd(i)
  482.     if(icop.eq.0)call
  483.      +     mulrat(itn(i),itd(i),in4,id4,itn(nt),itd(nt))
  484.     itb(nt)=k+1
  485.     do 10 j=ja,ite(i)
  486.     k=k+1
  487.     ivn(k)=ivn(j)
  488. 10    ivp(k)=ivp(j)
  489.     call nite(k)
  490. 15    continue
  491. 16    return
  492.     end
  493.  
  494.  
  495.     subroutine datim
  496. $include: 'picom.for'
  497.     integer*4 mult
  498.     call date(iy,mon,id)
  499.     call time(ih,min,is,if)
  500.     mult=60
  501.     itime=100*(is+mult*(min+mult*ih))+if
  502.     write(8,1)iy,mon,id,ih,min,is
  503. 1    format(' job run  ',i4,'/',i2,'/',i2,i9,':',i2,':',i2)
  504.     return
  505.     end
  506.  
  507.  
  508.     subroutine datime
  509. $include: 'picom.for'
  510.     integer*4 mult
  511.     call time(ih,min,is,if)
  512.     mult=60
  513.     itime=100*(is+mult*(min+mult*ih))+if
  514.     return
  515.     end
  516.  
  517.  
  518.     subroutine defh
  519. $include: 'picom.for'
  520.     character*1 ga(26),gc(26),gn(10)
  521.     data ga/'a','b','c','d','e','f','g','h','i','j','k','l','m'
  522.      ,       ,'n','o','p','q','r','s','t','u','v','w','x','y','z'/
  523.     data gc/'A','B','C','D','E','F','G','H','I','J','K','L','M'
  524.      ,       ,'N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/
  525.     data gn/'0','1','2','3','4','5','6','7','8','9'/
  526.     do 10 i=1,10
  527. 10    hn(i)=gn(i)
  528.     do 20 i=1,26
  529.     ha(i)=ga(i)
  530. 20    hc(i)=gc(i)
  531.     return
  532.     end
  533.  
  534.  
  535.     subroutine defsum(is,ittbv,ittev)
  536. $include: 'picom.for'
  537.     if(ittev.ne.nt .or. ittev+1.ne.ittbv)goto 1
  538.     call niltrm
  539.     ittev=nt
  540. 1    if(ittbv.gt.ittev)call bomb(3,'defsum1 ')
  541.     if(isb(is).eq.0)goto 20
  542.     j=ios(is)
  543.     minns=min0(minns,j)
  544.     if(ittbv.ge.isb(is) .and. ittev.le.ise(is))goto 8
  545.     if(j.eq.ns .and. ittbv.ge.isb(is))goto 8
  546.     if(ittbv.gt.ise(iosi(ns)))goto 4
  547.     call bomb(3,'defsum2 ')
  548. 4    ios(is)=ns
  549.     do 6 i=j,ns-1
  550.     k=iosi(i+1)
  551.     iosi(i)=k
  552. 6    ios(k)=i
  553.     iosi(ns)=is
  554. 8    isb(is)=ittbv
  555.     ise(is)=ittev
  556.     call pack
  557.     return
  558. 20    if(ittbv.le.ise(iosi(ns)))call bomb(3,'defsum3 ')
  559.     ns=ns+1
  560.     ios(is)=ns
  561.     iosi(ns)=is
  562.     isb(is)=ittbv
  563.     ise(is)=ittev
  564.     minns=min0(minns,ns)
  565.     call pack
  566.     return
  567.     end
  568.  
  569.  
  570.     subroutine delet
  571. $include: 'picom.for'
  572. c  delete sj4
  573.     call nxt(4,nchpl,' ',j1)
  574. 1    call numv(j1,nchpl,j2,j3,j4)
  575.     call delsum(j4)
  576.     if(j3.eq.nchpl)goto 10
  577.     call srchnb(j3+1,nchpl,j1)
  578.     if(j1.ne.-1)goto 1
  579. 10    return
  580.     end
  581.  
  582.  
  583.     subroutine delsum(is)
  584. $include: 'picom.for'
  585.     if(isb(is).eq.0)return
  586.     j=ios(is)
  587.     minns=min0(minns,j)
  588.     if(j.eq.0 .or. ns.eq.0)call bomb(3,'delsum1 ')
  589.     ns=ns-1
  590.     if(j.gt.ns)goto 10
  591.     do 2 i=j,ns
  592.     k=iosi(i+1)
  593.     iosi(i)=k
  594. 2    ios(k)=i
  595. 10    iosi(ns+1)=0
  596.     ios(is)=0
  597.     isb(is)=0
  598.     ise(is)=0
  599.     call pack
  600.     return
  601.     end
  602.  
  603.  
  604.     subroutine deriv(j10,ind,jans)
  605. c  d(sum j10)/d(var ind) is put in sum jans
  606. c  ider(1,nd) is dependent variable in derivative
  607. c  ider(2,nd) is independent variable in derivative
  608. c  ider(3,nd) is variable equal to derivative
  609. $include: 'picom.for'
  610.     integer*4 id1,ivpic
  611.     character*1 bl
  612.     common/i/if
  613.     data id1/1/
  614.     nt1=nt+1
  615.     if(isb(j10).lt.1 .or. ise(j10).gt.nt)call bomb(3,'deriv2  ')
  616.     do 90 it=isb(j10),ise(j10)
  617.     ia=itb(it)
  618.     ib=ite(it)
  619.     do 85 ic=ia,ib
  620.     ie=ivn(ic)
  621.     ivpic=ivp(ic)
  622.     if(ie.ne.ind)goto 40
  623.     call incnt
  624.     itb(nt)=itetop+1
  625.     call mulrat(itn(it),itd(it),ivpic,id1,itn(nt),itd(nt))
  626.     if=itetop
  627.     do 36 i=ia,ib
  628.     if(i.eq.ic)goto 32
  629.     if=if+1
  630.     ivn(if)=ivn(i)
  631.     ivp(if)=ivp(i)
  632.     goto 36
  633. 32    if(ivp(i).eq.1 .and. ia.ne.ib)goto 36
  634.     if=if+1
  635.     ivn(if)=ivn(i)
  636.     ivp(if)=ivp(i)-1
  637. 36    continue
  638.     call nite(if)
  639.     if(if.lt.itb(nt))call bomb(3,'deriv3  ')
  640.     goto 85
  641. 40    if(nd.eq.0)goto 85
  642.     i1=1
  643. 41    j1=idtab(i1,ie)
  644.     if(j1.eq.0)goto 85
  645.     k1=ider(2,j1)
  646.     level=1
  647.     if(k1.eq.ind)goto 50
  648.     i2=1
  649. 42    j2=idtab(i2,k1)
  650.     if(j2.eq.0)goto 77
  651.     k2=ider(2,j2)
  652.     level=2
  653.     if(k2.eq.ind)goto 50
  654.     i3=1
  655. 43    j3=idtab(i3,k2)
  656.     if(j3.eq.0)goto 76
  657.     k3=ider(2,j3)
  658.     level=3
  659.     if(k3.eq.ind)goto 50
  660.     i4=1
  661. 44    j4=idtab(i4,k3)
  662.     if(j4.eq.0)goto 75
  663.     k4=ider(2,j4)
  664.     level=4
  665.     if(k4.eq.ind)goto 50
  666.     i5=1
  667. 45    j5=idtab(i5,k4)
  668.     if(j5.eq.0)goto 74
  669.     k5=ider(2,j5)
  670.     level=5
  671.     if(k5.eq.ind)goto 50
  672.     l1=ipv(ie)
  673.     l2=ipv(k1)
  674.     l3=ipv(k2)
  675.     l4=ipv(k3)
  676.     l5=ipv(k4)
  677.     l6=ipv(k5)
  678.     bl=' '
  679.     write(8,49)(pv(i,ie),i=1,l1),bl,(pv(i,k1),i=1,l2),bl
  680.      +  ,(pv(i,k2),i=1,l3),bl,(pv(i,k3),i=1,l4),bl,(pv(i,k4),i=1,l5)
  681.      +  ,bl,(pv(i,k5),i=1,l6)
  682. 49    format(' chain too long: ',90a1)
  683.     call bomb(3,'deriv4  ')
  684. 50    call incnt
  685.     if=itetop
  686.     itb(nt)=if+1
  687.     call mulrat(itn(it),itd(it),ivpic,id1,itn(nt),itd(nt))
  688.     do 55 i=ia,ib
  689.     if(i.eq.ic)goto 52
  690.     if=if+1
  691.     ivn(if)=ivn(i)
  692.     ivp(if)=ivp(i)
  693.     goto 55
  694. 52    if(ivp(i).eq.1)goto 55
  695.     if=if+1
  696.     ivn(if)=ivn(i)
  697.     ivp(if)=ivp(i)-1
  698. 55    continue
  699.     call deriv2(j1)
  700.     if(level.eq.1)goto 77
  701.     call deriv2(j2)
  702.     if(level.eq.2)goto 76
  703.     call deriv2(j3)
  704.     if(level.eq.3)goto 75
  705.     call deriv2(j4)
  706.     if(level.eq.4)goto 74
  707.     call deriv2(j5)
  708.     i5=i5+1
  709.     if(i5.le.ndind)goto 45
  710. 74    i4=i4+1
  711.     if(i4.le.ndind)goto 44
  712. 75    i3=i3+1
  713.     if(i3.le.ndind)goto 43
  714. 76    i2=i2+1
  715.     if(i2.le.ndind)goto 42
  716. 77    i1=i1+1
  717.     if(i1.le.ndind)goto 41
  718. 85    continue
  719. 90    continue
  720.     call defsum(jans,nt1,nt)
  721.     call order(jans)
  722.     return
  723.     end
  724.  
  725.  
  726.     subroutine deriv2(j)
  727. $include: 'picom.for'
  728.     common/i/if
  729.     if=if+1
  730.     call nite(if)
  731.     ivn(if)=ider(3,j)
  732.     ivp(if)=1
  733.     return
  734.     end
  735.  
  736.  
  737.     subroutine dertab
  738. $include: 'picom.for'
  739.     idflg=0
  740.     do 302 i=1,ndind
  741.     do 302 j=1,nvar
  742. 302    idtab(i,j)=0
  743.     do 304 i=1,nd
  744.     idepv=ider(1,i)
  745.     do 303 j=1,ndind
  746.     if(idtab(j,idepv).eq.0)goto 304
  747. 303    continue
  748.     call bomb(4,'NDIND   ')
  749. 304    idtab(j,idepv)=i
  750.     return
  751.     end
  752.  
  753.  
  754.     subroutine dfine
  755. $include: 'picom.for'
  756.     call nxt(4,nchpl,' ',j1)
  757.     call nxtch(j1,nchpl,j2)
  758.     call numv(j2,nchpl,j4,j5,j6)
  759.     ittb6=nt+1
  760.     call srch(j5,nchpl,'=',j7)
  761.     if(j7.eq.-1)goto 410
  762.     call term(j7+1)
  763.     goto 412
  764. 410    call read(ib)
  765.     if(np.ge.2)call prntap
  766. 1    format(1x,127a1)
  767.     if(ib.eq.1)goto 412
  768.     call term(1)
  769.     goto 410
  770. 412    call defsum(1,ittb6,nt)
  771.     call order(1)
  772.     call outqu(1)
  773.     call swap(j6)
  774.     return
  775.     end
  776.  
  777.  
  778.     subroutine dump
  779. $include: 'picom.for'
  780.     character*1 bl
  781.     data bl/' '/
  782.     write(8,1)nd,nv,ns,np,nt,line,maxe,maxnt,maxite
  783. 1    format(/,' nd=',i3,' nv=',i3,' ns=',i3,' np=',i2,' nt=',i5,
  784.      +     ' line=',i3,' maxe=',i4,' maxnt=',i5,' maxite=',i6)
  785.     do 10 i=1,nv
  786.     j=ipv(i)
  787.     if(j.eq.nch)goto 3
  788.     jp=j+1
  789.     do 4 jj=jp,nch
  790. 4    pv(jj,i)=bl
  791. 3    write(8,5)i,(pv(k,i),k=1,nch),isb(i),ise(i),ios(i),iosi(i)
  792. 5    format(' i=',i5,'   pv=',15a1,'   isb=',i5,'   ise=',i5
  793.      +      ,'  ios=',i3,'  iosi=',i3)
  794. 10    continue
  795.     write(8,11)
  796. 11    format(' ')
  797.     write(8,13)
  798. 13    format('      i  itn  itd  itb  ite')
  799.     do 20 i=1,nt,4
  800.     j=i+1
  801.     k=i+2
  802.     l=i+3
  803.     write(8,21)i,itn(i),itd(i),itb(i),ite(i),j,itn(j),itd(j),itb(j),
  804.      +ite(j),k,itn(k),itd(k),itb(k),ite(k),l,itn(l),itd(l),itb(l),ite(l)
  805. 21    format(1x,4(i8,4i5))
  806. 20    continue
  807.     write(8,11)
  808.     do 22 i=1,nd
  809.     j=ider(3,i)
  810.     l=ipv(j)
  811. 22    write(8,24)ider(1,i),ider(2,i),(pv(m,j),m=1,l)
  812. 24    format(' the derivative of',i4,' wrt',i4,' is ', 15a1)
  813.     write(8,11)
  814.     write(8,40)(i,ivn(i),ivp(i),i=1,itetop)
  815. 40    format(10(1x,i5,2i3))
  816.     return
  817.     end
  818.  
  819.  
  820.     subroutine expnd(jb,jp,jans)
  821. $include: 'picom.for'
  822.     integer*4 kn,kd,ij,ik,jprvn,jprvd
  823.     ia=isb(jb)
  824.     if(ia.eq.0)call bomb(3,'expnd1  ')
  825.     if(itn(ia).ne.1 .or. itd(ia).ne.1)call bomb(3,'expnd2  ')
  826.     ic=itb(ia)
  827.     if(ivn(ic).ne.2 .or. ivp(ic).ne.0)call bomb(3,'expnd3  ')
  828.     ic=ia+1
  829.     id=ise(jb)
  830. 5    call copyt(ia,ia,1,1)
  831.     if(ic.le.id)goto 10
  832. 6    call defsum(jans,nt,nt)
  833.     return
  834. 10    ig=nt
  835.     ie=isb(jp)
  836.     kn=itn(ie)
  837.     kd=itd(ie)
  838.     if(kn.eq.0)goto 6
  839.     jprvb=ig
  840.     jprve=ig
  841.     jprvn=1
  842.     jprvd=1
  843.     do 20 kp=1,nb
  844.     if(kn.eq.0)goto 25
  845.     call mulrat(jprvn,jprvd,kn,kd*kp,jprvn,jprvd)
  846.     jtn(kp)=jprvn
  847.     jtd(kp)=jprvd
  848.     kn=kn-kd
  849.     jtb(kp)=nt+1
  850.     call mult(jprvb,jprve,ic,id)
  851.     jte(kp)=nt
  852.     if(jtb(kp).gt.nt)goto 25
  853.     jprvb=jtb(kp)
  854. 20    jprve=nt
  855.     kp=nb+1
  856. 25    kp=kp-1
  857.     do 40 i=1,kp
  858.     ih=jtb(i)
  859.     ii=jte(i)
  860.     ij=jtn(i)
  861.     ik=jtd(i)
  862.     do 30 il=ih,ii
  863. 30    call mulrat(itn(il),itd(il),ij,ik,itn(il),itd(il))
  864. 40    continue
  865.     call defsum(jans,ig,ii)
  866.     call order(jans)
  867.     return
  868.     end
  869.  
  870.  
  871.     subroutine extrct
  872. $include: 'picom.for'
  873.     integer*4 ii4
  874.     call nxt(2,nchpl,' ',j1)
  875.     call srch(j1,nchpl,'-',j6)
  876.     if(j6.ne.-1)goto 300
  877.     num=1
  878.     do 1 i=2,10
  879.     if(a(j1-1).eq.hn(i))num=i-1
  880. 1    continue
  881.     call numv(j1,nchpl,j3,j4,j5)
  882.     call nxt(j4,nchpl,'=',j6)
  883.     call nxt(j6,nchpl,'m',j7)
  884.     call numv(j7+1,nchpl,j8,j9,j10)
  885.     k=isb(j10)
  886.     if(k.eq.0)call bomb(3,'extrct1 ')
  887.     kk=itn(k)
  888.     if(kk.gt.0 .and. itd(k).eq.1)goto 15
  889. 12    write(8,14)itn(k),itd(k),(a(i),i=j8,j9)
  890. 14    format(' bad term number  ',2i12,5x, 15a1)
  891.     call bomb(3,'extrct2 ')
  892. 15    call nxt(j9+1,nchpl,'o',j11)
  893.     call numv(j11+2,nchpl,j12,j13,j14)
  894.     kl=ise(j14)+1-isb(j14)
  895.     if(kk.gt.kl)goto 12
  896.     km=isb(j14)+kk-1
  897.     nt1=nt+1
  898.     ktop=min0(km+num-1,ise(j14))
  899.     call copyt(km,ktop,1,1)
  900.     call defsum(j5,nt1,nt)
  901.     call outqu(j5)
  902.     return
  903. 300    call numv(j1,nchpl,j2,j3,j4)
  904.     call nxt(j3,j6,'=',j5)
  905.     call nxt(j5,j6,'t',j9)
  906.     call nxt(j9,j6,' ',j10)
  907.     call nxtnb(j10,j6,j11)
  908.     call srch(j11,j6,' ',j12)
  909.     if(j12.eq.-1)j12=j6
  910.     call getint(j11,j12-1,j13,ii4)
  911.     call srchnb(j6+1,nchpl,j8)
  912.     call srch(j8,nchpl,' ',j14)
  913.     call getint(j8,j14-1,j15,ii4)
  914.     call nxt(j14,nchpl,'o',j19)
  915.     if(a(j19+1).ne.'f' .and. a(j19+2).ne.' ')call bomb(3,'extrct5 ')
  916. 330    call numv(j19+2,nchpl,j16,j17,j18)
  917.     if(isb(j18).eq.0)call bomb(3,'extrct6 ')
  918.     if(j15.gt.ise(j18)-isb(j18)+1)j15=ise(j18)-isb(j18)+1
  919.     if(j15.lt.j13)call bomb(3,'extrct7 ')
  920.     nt1=nt+1
  921.     call copyt(isb(j18)+j13-1,isb(j18)+j15-1,1,1)
  922.     call defsum(j4,nt1,nt)
  923.     call outqu(j4)
  924.     return
  925.     end
  926.  
  927.  
  928.     subroutine fort(it,ib)
  929. $include: 'picom.for'
  930.     character*1 pa,aa
  931.     character*22 scr
  932.     common/pr1/pa(232),aa(22)
  933.     if(np.eq.0)return
  934.     if(itd(it).eq.1 .and. itn(it).eq.-1)then
  935.         ib=ib+1
  936.         a(ib)='-'
  937.         goto 42
  938.         endif
  939.     if(itd(it).eq.1 .and. itn(it).eq. 1)then
  940.         ib=ib+1
  941.         a(ib)='+'
  942.         goto 42
  943.         endif
  944.     write(scr,43)itn(it)
  945.     read(scr,9)(aa(kk),kk=1,20)
  946. 9    format(20a1)
  947.     j=0
  948.     do 10 i=1,20
  949.     if(aa(i).eq.' ')goto 10
  950.     if(j.eq.0 .and. aa(i).ne.'-')then
  951.         ib=ib+1
  952.         a(ib)='+'
  953.         endif
  954.     j=1
  955.     ib=ib+1
  956.     a(ib)=aa(i)
  957. 10    continue
  958.     ib=ib+1
  959.     a(ib)='.'
  960.     if(itd(it).eq.1)goto 42
  961.     ib=ib+1
  962.     a(ib)='/'
  963.     write(scr,43)itd(it)
  964.     read(scr,9)(aa(kk),kk=1,20)
  965.     do 15 i=1,20
  966.     if(aa(i).eq.' ')goto 15
  967.     ib=ib+1
  968.     a(ib)=aa(i)
  969. 15    continue
  970. 42    ic=itb(it)
  971.     id=ite(it)
  972.     do 60 iv=ic,id
  973.     ip=iabs(ivp(iv))
  974.     in=ivn(iv)
  975.     ie=ipv(in)
  976.     iq=0
  977.     if(ip.eq.1)goto 47
  978.     write(scr,43)ip
  979.     read(scr,9)(aa(kk),kk=1,20)
  980. 43    format(i20)
  981.     il=20
  982.     do 44 if=1,20
  983.     if(aa(if).ne.' ')goto 45
  984. 44    continue
  985. 45    iq=il+1-if+2
  986. 47    if(iq+ie+1.le.72)goto 52
  987.     do 48 j=75,1,-1
  988.     if(a(j).ne.' ')goto 49
  989. 48    continue
  990. 49    write(8,50)(a(i),i=1,j)
  991. 50    format(80a1)
  992.     do 51 i=1,80
  993. 51    a(i)=' '
  994.     a(6)='.'
  995.     ib=11
  996. 52    ib=ib+1
  997.     if(a(ib-1).ne.'+' .and. a(ib-1).ne.'-')a(ib)='*'
  998.     if(ivp(iv).lt.0)a(ib)='/'
  999.     do 53 i=1,ie
  1000.     ib=ib+1
  1001. 53    a(ib)=pv(i,in)
  1002.     if(iq.eq.0)goto 60
  1003.     ib=ib+1
  1004.     a(ib)='*'
  1005.     ib=ib+1
  1006.     a(ib)='*'
  1007.     do 54 i=if,il
  1008.     ib=ib+1
  1009. 54    a(ib)=aa(i)
  1010. 60    continue
  1011.     do 62 j=75,1,-1
  1012.     if(a(j).ne.' ')goto 63
  1013. 62    continue
  1014. 63    write(8,50)(a(i),i=1,j)
  1015.     return
  1016.     end
  1017.  
  1018.  
  1019.     subroutine fortrn
  1020. $include: 'picom.for'
  1021.     call nxt(1,nchpl,' ',j1)
  1022.     call numv(j1+1,nchpl,j2,j3,j4)
  1023.     if(np.eq.0 .or. isb(j4).eq.0)goto 10
  1024.     do 1 i=1,nchpl
  1025. 1    a(i)=' '
  1026.     ia=ipv(j4)
  1027.     ib=6
  1028.     do 2 i=1,ia
  1029.     ib=ib+1
  1030. 2    a(ib)=pv(i,j4)
  1031.     ib=ib+1
  1032.     a(ib)='='
  1033.     ic=isb(j4)
  1034.     id=ise(j4)
  1035.     j=0
  1036.     do 4 it=ic,id
  1037.     call fort(it,ib)
  1038.     do 3 i=1,nchpl
  1039. 3    a(i)=' '
  1040.     a(6)='.'
  1041.     ib=7
  1042.     j=j+1
  1043.     if(j.lt.10 .or. it.eq.id)goto 4
  1044.     a(6)=' '
  1045.     call copyc(pv(1,j4),a(7),ia)
  1046.     a(7+ia)='='
  1047.     call copyc(pv(1,j4),a(8+ia),ia)
  1048.     ib=7+ia+ia
  1049.     j=0
  1050. 4    continue
  1051. 10    return
  1052.     end
  1053.  
  1054.  
  1055.     subroutine gcd(inn,idd,ig)
  1056.     implicit integer*4 (i-n)
  1057.     ig=1
  1058.     in=iabs(inn)
  1059.     id=idd
  1060.     if(in.lt.id)goto 10
  1061. 5    in=mod(in,id)
  1062.     if(in.eq.0)goto 20
  1063. 10    if(in.eq.1)return
  1064.     id=mod(id,in)
  1065.     if(id.eq.0)goto 15
  1066.     if(id.eq.1)return
  1067.     goto 5
  1068. 15    ig=in
  1069.     return
  1070. 20    ig=id
  1071.     return
  1072.     end
  1073.  
  1074.  
  1075.     subroutine gcdout(in,id)
  1076.     implicit integer*4 (i-n)
  1077.     if(id.le.0)call bomb(3,'gcdout1 ')
  1078.     if(in.eq.1 .or. in.eq.-1 .or. id.eq.1)return
  1079.     if(in.eq.0)goto 1
  1080.     call gcd(in,id,ig)
  1081.     if(ig.eq.1)return
  1082.     in=in/ig
  1083.     id=id/ig
  1084.     return
  1085. 1    id=1
  1086.     return
  1087.     end
  1088.  
  1089.  
  1090.     subroutine getint(j1,j2,ii2,j3)
  1091. c  on return j3=integer value of a(j1)-a(j2)
  1092.     implicit integer*2 (i-n)
  1093.     integer*4 j3
  1094.     j3=0
  1095.     if(j2.lt.j1)call bomb(3,' getint1')
  1096.     do 1 i=j1,j2
  1097.     call intval(i,j)
  1098. 1    j3=10*j3+j
  1099.     ii2=j3
  1100.     return
  1101.     end
  1102.  
  1103.  
  1104.     subroutine getexp(j1,ie,k,if)
  1105. c  begin looking at a(j1), exponent put in ie, a(k) is next nonblank
  1106. c     (or k=-1 if only blanks are found)
  1107. $include: 'picom.for'
  1108.     ie=1
  1109.     call srchnb(j1,nchpl,k)
  1110.     if(k.eq.-1)return
  1111.     if(a(k).eq.'*' .and. a(k+1).ne.'*')return
  1112.     isgn=1
  1113.     if(a(k).eq.'*' .and. a(k+1).eq.'*')goto 35
  1114.     if(a(k).eq.'/')return
  1115.     call srchch(j1,nchpl,j2)
  1116.     if(j2.eq.k .and. if.ne.0)return
  1117.     if(j2.ne.k)goto 9
  1118. 32    call numv(j2,nchpl,j3,j4,j5)
  1119.     if(isb(j5).eq.0)call bomb(3,'getexp10')
  1120.     ie=itn(isb(j5))
  1121.     if(itd(isb(j5)).ne.1)call bomb(3,'getexp11')
  1122.     k=-1
  1123.     if(j4.ge.nchpl-1)return
  1124.     call srchnb(j4+1,nchpl,k)
  1125.     return
  1126. 9    if(a(k).eq.'+')goto 20
  1127.     if(a(k).eq.'-')goto 10
  1128.     if(a(k).eq.'^')goto 30
  1129.     j3=k
  1130. 1    me=0
  1131. 5    call intval(j3,j4)
  1132.     me=10*me+j4
  1133.     j3=j3+1
  1134.     if(j3.gt.nchpl)goto 70
  1135.     if(a(j3).eq.' ')goto 50
  1136.     if(a(j3).eq.'*')goto 60
  1137.     if(a(j3).eq.'/')goto 60
  1138.     goto 5
  1139. 10    isgn=-1
  1140. 20    if(k+1.gt.nchpl)call bomb(3,' getexp2')
  1141.     call nxtnb(k+1,nchpl,j3)
  1142.     if(j3.eq.-1)call bomb(3,' getexp4')
  1143.     goto 1
  1144. 30    call nxtnb(k+1,nchpl,j3)
  1145. 31    k=j3
  1146.     j2=j3
  1147.     call srchch(j3,nchpl,j4)
  1148.     if(j3.eq.j4)goto 32
  1149.     if(a(k).eq.'-')goto 10
  1150.     if(a(k).eq.'+')goto 20
  1151.     goto 1
  1152. 35    if(k.ge.nchpl-1)call bomb(3,' getexp6')
  1153.     call nxtnb(k+2,nchpl,j3)
  1154.     if(j3.eq.-1)call bomb(3,' getexp8')
  1155.     goto 31
  1156. 50    ie=me*isgn
  1157.     call srchnb(j3,nchpl,k)
  1158.     return
  1159. 60    k=j3
  1160.     if(j3.eq.nchpl)call bomb(3,' getexp9')
  1161.     ie=me*isgn
  1162.     return
  1163. 70    ie=me*isgn
  1164.     k=-1
  1165.     return
  1166.     end
  1167.  
  1168.  
  1169.     function icomp(ia,ib)
  1170. $include: 'picom.for'
  1171. c  if ia should precede ib icomp=1
  1172. c  if ib should precede ia icomp=-1
  1173. c  if terms are the same except possibly coefficients, icomp=0
  1174.     ja=itb(ia)
  1175.     jb=itb(ib)
  1176.     ka=ite(ia)
  1177.     kb=ite(ib)
  1178. 1    if(ivn(ja)-ivn(jb))2,8,9
  1179. 2    if(ivp(ja))3,10,4
  1180. 4    icomp=-1
  1181.     return
  1182. 5    if(jb.gt.kb)goto 7
  1183. 3    icomp=1
  1184.     return
  1185. 6    ja=ja+1
  1186. 11    jb=jb+1
  1187.     if(ja.gt.ka)goto 5
  1188.     if(jb.gt.kb)goto 4
  1189.     goto 1
  1190. 7    icomp=0
  1191.     return
  1192. 8    if(ivp(ja)-ivp(jb))3,6,4
  1193. 9    if(ivp(jb))4,11,3
  1194. 10    ja=ja+1
  1195.     if(ja.gt.ka)goto 3
  1196.     goto 1
  1197.     end
  1198.  
  1199.  
  1200.     subroutine iftest(igoflg)
  1201. $include: 'picom.for'
  1202.     integer*4 i5,i8
  1203.     character*1 aj1,aj2
  1204.     igoflg=0
  1205.     call nxt(2,nchpl,' ',j1)
  1206.     call nxt(j1,nchpl,'.',j)
  1207.     call numv(j1,j-1,j3,j4,j5)
  1208.     aj1=a(j+1)
  1209.     aj2=a(j+2)
  1210.     if(aj1.eq.'n')goto 20
  1211.     if(aj1.eq.'l' .or. aj1.eq.'g')goto 22
  1212.     if(aj1.ne.'e')call bomb(3,' iftest1')
  1213.     if(aj2.ne.'q')call bomb(3,' iftest2')
  1214. 4    if(a(j+3).ne.'.')call bomb(3,' iftest3')
  1215.     call numv(j+4,nchpl,j6,j7,j8)
  1216.     if(j5.ne.j8)goto 10
  1217.     if(aj1.eq.'n')return
  1218. 5    call nxtch(j7+1,nchpl,j9)
  1219.     do 6 j=j9,nchpl
  1220.     a(j+1-j9)=a(j)
  1221. 6    a(j)=' '
  1222.     igoflg=1
  1223.     call copyc(a,ap,nchpl)
  1224.     return
  1225. 10    k1=itb(isb(j5))
  1226.     k2=ite(ise(j5))
  1227.     k3=itb(isb(j8))
  1228.     k4=ite(ise(j8))
  1229.     i5=itn(isb(j5))*itd(isb(j8))
  1230.     i8=itd(isb(j5))*itn(isb(j8))
  1231.     if(aj1.eq.'l')goto 23
  1232.     if(aj1.eq.'g')goto 25
  1233.     if(k2-k1.ne.k4-k3)goto 30
  1234.     k=k3
  1235.     do 11 i=k1,k2
  1236.     if(ivn(i).ne.ivn(k))goto 30
  1237.     if(ivp(i).ne.ivp(k))goto 30
  1238. 11    k=k+1
  1239.     k1=isb(j5)
  1240.     k2=ise(j5)
  1241.     k3=isb(j8)
  1242.     k4=ise(j8)
  1243.     if(k2-k1.ne.k4-k3)goto 30
  1244.     k=k3
  1245.     do 12 i=k1,k2
  1246.     if(itn(i).ne.itn(k))goto 30
  1247.     if(itd(i).ne.itd(k))goto 30
  1248. 12    k=k+1
  1249.     if(aj1.eq.'n')return
  1250.     goto 5
  1251. 20    if(a(j+2).ne.'e')call bomb(3,' iftest4')
  1252.     goto 4
  1253. 22    if(aj2.eq.'e' .or. aj2.eq.'t')goto 4
  1254.     call bomb(5,' iftest5')
  1255. 23    if(aj2.eq.'e')goto 24
  1256.     if(i5.lt.i8)goto 5
  1257.     return
  1258. 24    if(i5.le.i8)goto 5
  1259.     return
  1260. 25    if(aj2.eq.'e')goto 26
  1261.     if(i5.gt.i8)goto 5
  1262.     return
  1263. 26    if(i5.ge.i8)goto 5
  1264.     return
  1265. 30    if(aj1.eq.'n')goto 5
  1266.     return
  1267.     end
  1268.  
  1269.  
  1270.     subroutine incnt
  1271. $include: 'picom.for'
  1272.     nt=nt+1
  1273.     maxnt=max0(nt,maxnt)
  1274.     if(maxnt.gt.nterm)call bomb(4,'NTERM   ')
  1275.     return
  1276.     end
  1277.  
  1278.  
  1279.     subroutine intval(j1,i)
  1280. c  on return i=integer value of a(j1)
  1281. $include: 'picom.for'
  1282.     do 1 j=1,10
  1283.     if(a(j1).eq.hn(j))goto 2
  1284. 1    continue
  1285.     i=-1
  1286.     return
  1287. 2    i=j-1
  1288.     return
  1289.     end
  1290.  
  1291.  
  1292.     subroutine monitr
  1293. $include: 'picom.for'
  1294.     call datime
  1295.     tprev=itime/100.
  1296.     call readin
  1297.     ntp=0
  1298. 10    if(ha(1).ne.'a')call bomb(3,'monitr1 ')
  1299.     if(isb(2).ne.0)call bomb(3,'useorder')
  1300.     if(ha(1).eq.'a')goto 11
  1301.     if(ns.lt.2)goto 14
  1302.     do 20 ms=2,ns
  1303.     is=iosi(ms-1)
  1304.     js=iosi(ms)
  1305.     if(is.eq.0 .or. js.eq.0)call bomb(3,'    mon1')
  1306.     if(ios(is).ne.ms-1 .or. ios(js).ne.ms)call bomb(3,'    mon2')
  1307.     if(isb(is).eq.0 .or. isb(js).eq.0)call bomb(3,'    mon3')
  1308.     ia=ise(is)
  1309.     ib=isb(js)
  1310.     if(ite(ia)+1.ne.itb(ib))call bomb(3,'    mon4')
  1311.     do 16 it=ib,ise(js)
  1312.     if(itb(it).gt.ite(it))call bomb(3,'    mon5')
  1313.     if(itd(it).le.0)call bomb(3,'   mon22')
  1314.     if(itn(it).ne.0)goto 15
  1315.     if(itb(it).ne.ite(it))call bomb(3,'    mon6')
  1316.     if(ivn(itb(it)).ne.2)call bomb(3,'    mon7')
  1317.     if(itb(it).gt.ite(it))call bomb(3,'    mon8')
  1318. 15    do 17 i=itb(it),ite(it)
  1319.     if(ivn(i).eq.1)call bomb(3,'    mon9')
  1320.     if(isb(ivn(i)).ne.0)call bomb(3,'   mon10')
  1321.     if(i.eq.itb(it))goto 17
  1322.     if(ivp(i).eq.0)call bomb(3,'   mon11')
  1323.     if(ivn(i-1).ge.ivn(i))call bomb(3,'   mon12')
  1324. 17    continue
  1325. 16    continue
  1326. 20    continue
  1327. 14    if(ns.eq.0)call bomb(3,'   mon13')
  1328.     ic=iosi(ns)
  1329.     if(ise(ic).ne.nt)call bomb(3,'   mon14')
  1330.     if(ite(ise(ic)).ne.itetop)call bomb(3,'   mon15')
  1331. 11    call read(ib)
  1332.     if(ib.eq.1 .and. np.ge.1)call prntap
  1333. 1    format(1x,127a1)
  1334.     if(ib.eq.1)goto 10
  1335.     call datime
  1336.     ttot=(itime-itbeg)/100.
  1337.     t=ttot-tprev
  1338.     tprev=ttot
  1339.     ntmntp=nt-ntp
  1340.     if(np.ge.3)write(8,3)t,ttot,ntmntp,nt,itetop,maxnt,maxite
  1341. 3    format(10x,'t=',f7.2,' total t=',f9.1,' ntmntp=',i5,' nt=',i5,
  1342.      +       ' itetop=',i6,' maxnt=',i5,' maxite=',i6)
  1343.     if(np.ge.2)call prntap
  1344. 5    ntp=nt
  1345.     call copyc(a,ap,nchpl)
  1346.     if(a(1).eq.'c')goto 300
  1347.     if(a(1).eq.'d')goto 400
  1348.     if(a(1).eq.'e')goto 500
  1349.     if(a(1).eq.'f')goto 600
  1350.     if(a(1).eq.'g')goto 700
  1351.     if(a(1).eq.'i')goto 900
  1352.     if(a(1).eq.'l')goto 10
  1353.     if(a(1).eq.'m')goto 1300
  1354.     if(a(1).eq.'n')goto 1400
  1355.     if(a(1).eq.'o')goto 1500
  1356.     if(a(1).eq.'p')goto 1600
  1357.     if(a(1).eq.'r')goto 1800
  1358.     if(a(1).eq.'s')goto 1900
  1359.     if(a(1).eq.'t')goto 2000
  1360.     if(a(1).eq.'z')goto 2600
  1361.     call bomb(3,'bad cmnd')
  1362. 300    if(a(2).eq.'o')goto 310
  1363.     if(a(2).ne.'a')call bomb(3,' monitr2')
  1364.     call call
  1365.     goto 10
  1366. 310    if(a(3).eq.'m')goto 320
  1367.     if(a(3).ne.'u')call bomb(3,' monitr3')
  1368.     call count
  1369.     goto 10
  1370. 320    call comput
  1371.     call pack
  1372.     goto 10
  1373. 400    if(a(2).eq.'u')goto 470
  1374.     if(a(2).eq.'i')goto 320
  1375.     if(a(2).ne.'e')call bomb(3,' monitr4')
  1376.     if(a(3).eq.'l')goto 450
  1377.     if(a(3).eq.'f')goto 440
  1378.     call bomb(3,' monitr5')
  1379. 440    call dfine
  1380.     goto 10
  1381. 450    call delet
  1382.     call pack
  1383.     goto 10
  1384. 470    call dump
  1385.     goto 10
  1386. 500    if(a(2).eq.'n')return
  1387.     if(a(2).ne.'x')call bomb(3,' monitr6')
  1388.     call extrct
  1389.     goto 10
  1390. 600    if(a(2).eq.'l')goto 610
  1391.     call fortrn
  1392.     goto 10
  1393. 610    call readin
  1394.     goto 10
  1395. 700    call nxt(1,nchpl,' ',j9)
  1396.     call numl(j9,nchpl,j10,j11,j12)
  1397.     if(linlab(j12).eq.-1)write(8,710)
  1398.     if(linlab(j12).eq.-1)call bomb(3,' monitr7')
  1399. 710    format(' label not found')
  1400.     line=linlab(j12)-1
  1401.     goto 10
  1402. 900    call iftest(igoflg)
  1403.     if(igoflg.eq.1)goto 5
  1404.     goto 10
  1405. 1300    call srch(1,nchpl,'=',j)
  1406.     call srchnb(j+1,nchpl,k)
  1407.     call getexp(k,maxe,j,0)
  1408.     goto 10
  1409. 1400    call nxt(3,nchpl,'=',j1)
  1410.     call getexp(j1+1,j3,j2,0)
  1411.     if(a(2).eq.'d')goto 1405
  1412.     if(a(2).eq.'p')goto 1410
  1413.     if(a(2).eq.'b')goto 1415
  1414.     call bomb(3,' monitr8')
  1415. 1405    if(a(3).eq.'u')goto 1406
  1416.     call bomb(3,' monitr9')
  1417. 1406    ndump=j3
  1418.     goto 10
  1419. 1410    np=j3
  1420.     goto 10
  1421. 1415    nb=j3
  1422.     goto 10
  1423. 1500    if(nv.ne.1)call bomb(3,'putfirst')
  1424.     call ordset
  1425.     goto 10
  1426. 1600    if(a(3).eq.'i')goto 1620
  1427.     if(a(3).ne.'o')call bomb(3,'monitr10')
  1428. 1605    call read(ib)
  1429.     if(np.ge.2)call prntap
  1430.     if(a(1).ne.'r' .or. a(2).ne.'e' .or. a(3).ne.'t')goto 1605
  1431.     goto 10
  1432. 1620    call nxt(4,nchpl,' ',j1)
  1433.     call numv(j1,nchpl,j2,j3,j4)
  1434.     if(np.eq.1)write(8,1)
  1435.     if(np.eq.1)call prntap
  1436.     call print(j4)
  1437.     goto 10
  1438. 1800    if(a(2).ne.'e' .or. a(3).ne.'t')call bomb(3,' return?')
  1439.     line=linret
  1440.     goto 10
  1441. 1900    call sub
  1442.     call pack
  1443.     goto 10
  1444. 2000    if(a(2).eq.'a')goto 2010
  1445.     call thderv
  1446.     goto 10
  1447. 2010    call taylor
  1448.     goto 10
  1449. 2600    call numv(2,nchpl,j1,j2,j3)
  1450.     call zsub(j3)
  1451.     call pack
  1452.     goto 10
  1453.     end
  1454.  
  1455.  
  1456.     subroutine mulrat(n1,i1,n2,i2,n3,i3)
  1457.     implicit integer*4 (i-n)
  1458.     n3=n1*n2
  1459.     i3=i1*i2
  1460.     call gcdout(n3,i3)
  1461.     return
  1462.     end
  1463.  
  1464.  
  1465.     subroutine mult(ia,ib,ic,id)
  1466. c  multiply terms ia-ib by terms ic-id
  1467. c  result goes to nt+1,.....
  1468. $include: 'picom.for'
  1469.     common/nwt/iitb,iite,jitb,jite,ntb,if
  1470.     integer*4 inp,idp
  1471.     if(ia.eq.0 .or. ic.eq.0)call bomb(3,'   mult1')
  1472.     ntb=nt+1
  1473.     if=0
  1474.     do 20 i=ia,ib
  1475.     iitb=itb(i)
  1476.     iite=ite(i)
  1477.     do 10 j=ic,id
  1478.     jitb=itb(j)
  1479.     jite=ite(j)
  1480.     if(ivn(iitb)+ivn(jitb).eq.4 .and. ivp(iitb)+ivp(jitb).gt.maxe)
  1481.      +                                           goto 20
  1482.     call mulrat(itn(i),itd(i),itn(j),itd(j),inp,idp)
  1483.     call sub2(inp,idp)
  1484. 10    continue
  1485. 20    continue
  1486.     return
  1487.     end
  1488.  
  1489.  
  1490.     subroutine multn(i,j,k,l)
  1491. $include: 'picom.for'
  1492.     nt1=nt+1
  1493.     call mult(i,j,k,l)
  1494.     if(nt.ge.nt1)return
  1495.     call niltrm
  1496.     return
  1497.     end
  1498.  
  1499.  
  1500.     subroutine niltrm
  1501. $include: 'picom.for'
  1502.     k=itetop+1
  1503.     call incnt
  1504.     itn(nt)=0
  1505.     itd(nt)=1
  1506.     itb(nt)=k
  1507.     call nite(k)
  1508.     ivn(k)=2
  1509.     ivp(k)=0
  1510.     return
  1511.     end
  1512.  
  1513.  
  1514.     subroutine nite(j)
  1515. $include: 'picom.for'
  1516.     itetop=j
  1517.     ite(nt)=j
  1518.     maxite=max0(maxite,j)
  1519.     if(maxite.gt.ntot)call bomb(4,'NTOT    ')
  1520.     return
  1521.     end
  1522.  
  1523.  
  1524.  
  1525.     subroutine numl(j1,j2,j3,j4,j5)
  1526. c  search j1-j2  find name in j3-j4, j5 is index of label
  1527. $include: 'picom.for'
  1528.     character*1 c
  1529.     if(j2.lt.j1)call bomb(3,'    numl')
  1530.     call nxtch(j1,j2,j3)
  1531.     do 10 j=j3,j2
  1532.     c=a(j)
  1533.     do 5 k=1,26
  1534.     if(c.eq.ha(k) .or. c.eq.hc(k))goto 10
  1535. 5    continue
  1536.     do 6 k=1,10
  1537.     if(c.eq.hn(k))goto 10
  1538. 6    continue
  1539.     if(c.eq.' ' .or. c.eq.'+' .or. c.eq.'-')goto 12
  1540.     if(c.eq.'*' .or. c.eq.'/' .or. c.eq.'=' .or. c.eq.'^')goto 12
  1541. 10    continue
  1542.     j4=j2
  1543.     goto 14
  1544. 12    j4=j-1
  1545. 14    l=j4-j3+1
  1546.     do 18 i=1,labels
  1547.     if(ipl(i).ne.l)goto 18
  1548.     do 16 j=1,l
  1549.     if(pl(j,i).ne.a(j+j3-1))goto 18
  1550. 16    continue
  1551.     j5=i
  1552.     return
  1553. 18    continue
  1554.     labels=labels+1
  1555.     if(labels.gt.nlab)call bomb(4,'NLAB    ')
  1556.     j5=labels
  1557.     ipl(j5)=l
  1558.     call copyc(a(j3),pl(1,j5),l)
  1559.     return
  1560.     end
  1561.  
  1562.  
  1563.     subroutine numv(j1,j2,j3,j4,j5)
  1564. c  search j1-j2  find name in j3-j4, j5 is index of variable
  1565. $include: 'picom.for'
  1566.     character*1 c
  1567.     if(j2.lt.j1)call bomb(3,'    numv')
  1568.     call nxtch(j1,j2,j3)
  1569.     do 10 j=j3,j2
  1570.     c=a(j)
  1571.     do 5 k=1,26
  1572.     if(c.eq.ha(k) .or. c.eq.hc(k))goto 10
  1573. 5    continue
  1574.     do 6 k=1,10
  1575.     if(c.eq.hn(k))goto 10
  1576. 6    continue
  1577.     if(c.eq.' ' .or. c.eq.'+' .or. c.eq.'-')goto 12
  1578.     if(c.eq.'*' .or. c.eq.'/' .or. c.eq.'=' .or. c.eq.'^')goto 12
  1579. 10    continue
  1580.     j4=j2
  1581.     goto 14
  1582. 12    j4=j-1
  1583. 14    l=j4-j3+1
  1584.     do 28 ii=ns,1,-1
  1585.     i=iosi(ii)
  1586.     if(ipv(i).ne.l)goto 28
  1587.     do 26 j=1,l
  1588.     if(pv(j,i).ne.a(j+j3-1))goto 28
  1589. 26    continue
  1590.     j5=i
  1591.     return
  1592. 28    continue
  1593.     do 18 i=1,nv
  1594.     if(isb(i).ne.0)goto 18
  1595.     if(ipv(i).ne.l)goto 18
  1596.     do 16 j=1,l
  1597.     if(pv(j,i).ne.a(j+j3-1))goto 18
  1598. 16    continue
  1599.     j5=i
  1600.     return
  1601. 18    continue
  1602.     nv=nv+1
  1603.     j5=nv
  1604.     if(nv.gt.nvar)call bomb(4,'NVAR    ')
  1605.     ipv(nv)=l
  1606.     call copyc(a(j3),pv(1,nv),l)
  1607.     return
  1608.     end
  1609.  
  1610.  
  1611.     subroutine nxt(j1,j2,h,j3)
  1612.     implicit integer*2 (i-n)
  1613.     character*1 h
  1614.     call srch(j1,j2,h,j3)
  1615.     if(j3.eq.-1)call bomb(3,'     nxt')
  1616.     end
  1617.  
  1618.  
  1619.     subroutine nxtch(j1,j2,j3)
  1620.     implicit integer*2 (i-n)
  1621.     call srchch(j1,j2,j3)
  1622.     if(j3.eq.-1)call bomb(3,'   nxtch')
  1623.     return
  1624.     end
  1625.  
  1626.  
  1627.     subroutine nxtnb(j1,j2,j3)
  1628.     implicit integer*2 (i-n)
  1629.     call srchnb(j1,j2,j3)
  1630.     if(j3.eq.-1)call bomb(3,'   nxtnb')
  1631.     end
  1632.  
  1633.  
  1634.     subroutine nxtnnu(j1,j2,j3)
  1635. $include: 'picom.for'
  1636.     do 10 i=j1,j2
  1637.     do 5 j=1,10
  1638.     if(a(i).eq.hn(j))goto 10
  1639. 5    continue
  1640.     j3=i
  1641.     return
  1642. 10    continue
  1643.     j3=-1
  1644.     return
  1645.     end
  1646.  
  1647.  
  1648.     subroutine order(is)
  1649. $include: 'picom.for'
  1650.     ina=isb(is)
  1651.     inb=ise(is)
  1652.     if(inb.lt.ina .or. ina.eq.0)call bomb(3,'  order1')
  1653.     it=ina
  1654. 2    ia=itb(it)
  1655.     if(itn(it).eq.0)goto 41
  1656.     ib=ite(it)
  1657. 3    if(ib.eq.ia)goto 25
  1658.     i=ia
  1659. 4    if(ivp(i).eq.0)goto 20
  1660.     i=i+1
  1661.     if(i.le.ib)goto 4
  1662.     if(ia.eq.ib)goto 25
  1663. 23    i=ia+1
  1664. 5    if(ivn(i)-ivn(i-1))10,7,6
  1665. 6    i=i+1
  1666.     if(i.le.ib)goto 5
  1667.     goto 40
  1668. 7    ivp(i-1)=ivp(i-1)+ivp(i)
  1669.     if(i.eq.ib)goto 9
  1670.     ib=ib-1
  1671.     do 8 j=i,ib
  1672.     ivn(j)=ivn(j+1)
  1673. 8    ivp(j)=ivp(j+1)
  1674.     if(ivp(i-1).ne.0)goto 5
  1675.     goto 3
  1676. 9    ib=ib-1
  1677.     if(ivp(ib).ne.0)goto 40
  1678.     goto 3
  1679. 10    if(i.eq.ia+1)goto 19
  1680.     ic=i-2
  1681. 11    if(ivn(i)-ivn(ic))12,13,16
  1682. 12    if(ic.eq.ia)goto 18
  1683.     ic=ic-1
  1684.     goto 11
  1685. 13    ivp(ic)=ivp(ic)+ivp(i)
  1686.     if(i.eq.ib)goto 15
  1687.     ib=ib-1
  1688.     do 14 j=i,ib
  1689.     ivn(j)=ivn(j+1)
  1690. 14    ivp(j)=ivp(j+1)
  1691.     if(ivp(ic).eq.0)goto 3
  1692.     goto 5
  1693. 15    ib=ib-1
  1694.     if(ivp(ic).ne.0)goto 40
  1695.     goto 3
  1696. 16    ip=ivp(i)
  1697.     in=ivn(i)
  1698.     do 17 j=i,ic+2,-1
  1699.     ivn(j)=ivn(j-1)
  1700. 17    ivp(j)=ivp(j-1)
  1701.     ivn(ic+1)=in
  1702.     ivp(ic+1)=ip
  1703.     goto 6
  1704. 18    ic=ia-1
  1705.     goto 16
  1706. 19    ip=ivp(i)
  1707.     in=ivn(i)
  1708.     ivp(i)=ivp(ia)
  1709.     ivn(i)=ivn(ia)
  1710.     ivp(ia)=ip
  1711.     ivn(ia)=in
  1712.     goto 6
  1713. 20    if(i.eq.ib)goto 22
  1714.     ib=ib-1
  1715.     do 21 j=i,ib
  1716.     ivn(j)=ivn(j+1)
  1717. 21    ivp(j)=ivp(j+1)
  1718.     goto 4
  1719. 22    ib=ib-1
  1720.     if(ib.gt.ia)goto 23
  1721.     ib=ia
  1722. 25    if(ivp(ia).eq.0)ivn(ia)=2
  1723. 40    ite(it)=ib
  1724.     if(ivn(ia).ne.2 .or. ivp(ia).le.maxe)goto 52
  1725. 41    call bkup(it,inb)
  1726.     goto 53
  1727. 52    it=it+1
  1728. 53    if(it.le.inb)goto 2
  1729. c  sort terms ina-inb
  1730.     if(ina.lt.inb)goto 70
  1731.     if(ina.eq.inb)goto 82
  1732. 51    inb=ina
  1733.     itn(ina)=0
  1734.     itd(ina)=1
  1735.     k=itb(ina)
  1736.     ite(ina)=k
  1737.     ivn(k)=2
  1738.     ivp(k)=0
  1739.     goto 82
  1740. 70    it=ina
  1741. 68    it=it+1
  1742. 69    if(it.gt.inb)goto 80
  1743.     if(itn(it).eq.0)goto 75
  1744.     if(icomp(ina,it))71,73,74
  1745. 71    call tzap(it,ina)
  1746.     goto 68
  1747. 73    call addrat(itn(ina),itd(ina),itn(it),itd(it),itn(ina),itd(ina))
  1748. 75    call bkup(it,inb)
  1749.     goto 69
  1750. 74    if(icomp(it-1,it))62,61,68
  1751. 61    call addrat(itn(it-1),itd(it-1),itn(it),itd(it),
  1752.      ,     itn(it-1),itd(it-1))
  1753.     call bkup(it,inb)
  1754.     goto 69
  1755. 62    la=ina
  1756.     lb=it-1
  1757. 63    if(lb-la.eq.1)goto 64
  1758.     lc=(la+lb)/2
  1759.     if(icomp(lc,it))65,66,67
  1760. 64    call tzap(it,lb)
  1761.     goto 68
  1762. 65    lb=lc
  1763.     goto 63
  1764. 66    call addrat(itn(lc),itd(lc),itn(it),itd(it),itn(lc),itd(lc))
  1765.     call bkup(it,inb)
  1766.     goto 69
  1767. 67    la=lc
  1768.     goto 63
  1769. 80    if(ina.gt.inb)goto 51
  1770. 82    call defsum(is,ina,inb)
  1771.     return
  1772.     end
  1773.  
  1774.  
  1775.     subroutine ordset
  1776. $include: 'picom.for'
  1777.     call nxt(2,nchpl,' ',j1)
  1778.     call srchnb(j1,nchpl,j2)
  1779.     if(j2.eq.-1)goto 20
  1780.     ntsav=nt
  1781.     itesav=itetop
  1782.     call term(j1)
  1783.     nt=ntsav
  1784.     itetop=itesav
  1785.     return
  1786. 20    if=1
  1787.     do 10 il=1,300
  1788.     call read(ib)
  1789.     if(np.ge.2)call prntap
  1790. 1    format(1x,127a1)
  1791.     if(ib.eq.1)return
  1792.     is=1
  1793. 2    call srchch(is,nchpl,j1)
  1794.     if(j1.eq.-1)goto 10
  1795.     call numv(j1,nchpl,j2,j3,j4)
  1796.     if=if+1
  1797.     if(j4.lt.if)call bomb(3,'ordset1 ')
  1798.     if(j4.eq.if)goto 8
  1799.     do 5 i=nv,if,-1
  1800.     ipv(i)=ipv(i-1)
  1801.     do 4 ic=1,nch
  1802. 4    pv(ic,i)=pv(ic,i-1)
  1803. 5    continue
  1804.     call copyc(a(j2),pv(1,if),j3-j2+1)
  1805.     ipv(if)=j3-j2+1
  1806. 8    is=j3+1
  1807.     if(is.le.nchpl)goto 2
  1808. 10    continue
  1809.     return
  1810.     end
  1811.  
  1812.  
  1813.     subroutine outqu(j5)
  1814. $include: 'picom.for'
  1815. 1    do 6 it=isb(j5),ise(j5)
  1816.     do 5 iv=itb(it),ite(it)
  1817.     if(isb(ivn(iv)).ne.0)goto 9
  1818. 5    continue
  1819. 6    continue
  1820.     return
  1821. 9    jv=ivn(iv)
  1822.     call subs(jv,jv,j5)
  1823.     goto 1
  1824.     end
  1825.  
  1826.  
  1827.     subroutine pack
  1828. $include: 'picom.for'
  1829.     if(minns.le.ns)goto 9
  1830.     is=iosi(ns)
  1831.     nt=ise(is)
  1832.     itetop=ite(nt)
  1833.     minns=ns+1
  1834.     return
  1835. 9    if(minns.le.0)call bomb(3,'   pack4')
  1836.     js=minns
  1837.     kt=0
  1838.     kv=0
  1839.     if(js.eq.1)goto 13
  1840.     kt=ise(iosi(js-1))
  1841.     kv=ite(kt)
  1842. 13    do 30 ms=js,ns
  1843.     is=iosi(ms)
  1844.     if(is.le.0)call bomb(3,'   pack6')
  1845.     ia=isb(is)
  1846.     ib=ise(is)
  1847.     lv=0
  1848.     lt=0
  1849.     do 14 it=ia,ib
  1850.     if(itn(it).eq.0)goto 14
  1851.     lt=lt+1
  1852.     if(lt.gt.ntermt)call bomb(4,'NTERMT  ')
  1853.     jtn(lt)=itn(it)
  1854.     jtd(lt)=itd(it)
  1855.     jtb(lt)=lv+1
  1856.     do 12 i=itb(it),ite(it)
  1857.     if(ivp(i).ne.0)goto 31
  1858.     do 32 j=itb(it),ite(it)
  1859.     if(ivp(j).ne.0)goto 12
  1860. 32    continue
  1861.     lv=lv+1
  1862.     jvn(lv)=2
  1863.     jvp(lv)=0
  1864.     goto 33
  1865. 31    lv=lv+1
  1866.     jvn(lv)=ivn(i)
  1867.     jvp(lv)=ivp(i)
  1868. 12    continue
  1869. 33    jte(lt)=lv
  1870.     if(lv.gt.ntott)call bomb(4,'NTOTT   ')
  1871. 14    continue
  1872.     if(lv.ne.0)goto 15
  1873.     kt=kt+1
  1874.     isb(is)=kt
  1875.     ise(is)=kt
  1876.     kv=kv+1
  1877.     itn(kt)=0
  1878.     itd(kt)=1
  1879.     itb(kt)=kv
  1880.     ite(kt)=kv
  1881.     ivn(kv)=2
  1882.     ivp(kv)=0
  1883.     goto 30
  1884. 15    isb(is)=kt+1
  1885.     do 18 it=1,lt
  1886.     kt=kt+1
  1887.     itn(kt)=jtn(it)
  1888.     itd(kt)=jtd(it)
  1889.     itb(kt)=kv+1
  1890.     do 16 i=jtb(it),jte(it)
  1891.     kv=kv+1
  1892.     ivn(kv)=jvn(i)
  1893. 16    ivp(kv)=jvp(i)
  1894. 18    ite(kt)=kv
  1895.     ise(is)=kt
  1896. 30    continue
  1897.     nt=kt
  1898.     itetop=kv
  1899.     minns=ns+1
  1900.     return
  1901.     end
  1902.  
  1903.  
  1904.     subroutine print(ii)
  1905. $include: 'picom.for'
  1906.     character*127 scr
  1907.     character*1 pa,aa,bl,sm,cf
  1908.     common/pr1/pa(232),aa(22)
  1909.     data bl,sm,cf/' ','-','^'/
  1910.     if(np.eq.0)return
  1911.     ia=isb(ii)
  1912.     if(ia.eq.0)write(8,2)
  1913.     if(ia.eq.0)return
  1914. 2    format(' there is nothing to print')
  1915.     ib=ise(ii)
  1916.     nst=23
  1917.     do 40 it=ia,ib
  1918.     n=nst
  1919.     do 4 i=1,232
  1920. 4    pa(i)=' '
  1921.     ka=itb(it)
  1922.     kb=ite(it)
  1923.     do 30 k=ka,kb
  1924.     n=n+1
  1925.     pa(n)=bl
  1926.     if=ivn(k)
  1927.     jp=ivp(k)
  1928.     ifl=ipv(if)
  1929.     do 20 j=1,ifl
  1930.     n=n+1
  1931. 20    pa(n)=pv(j,if)
  1932.     if(jp.eq.1)goto 30
  1933.     n=n+1
  1934.     if(jp.ge.0)pa(n)=cf
  1935.     if(jp.lt.0)pa(n)=sm
  1936.     iap=iabs(jp)
  1937.     if(iap.le.99)goto 21
  1938.     pa(n+1)='*'
  1939.     pa(n+2)='*'
  1940.     n=n+2
  1941.     goto 30
  1942. 21    n=n+1
  1943.     j=mod(iap,10)
  1944.     pa(n)=hn(j+1)
  1945.     if(iap.le.9)goto 30
  1946.     j=iap/10
  1947.     n=n+1
  1948.     pa(n)=pa(n-1)
  1949.     pa(n-1)=hn(j+1)
  1950. 30    if(n.gt.232)call bomb(3,'  print3')
  1951.     write(scr,31)itn(it)
  1952.     read(scr,36)(pa(kk),kk=1,11)
  1953. 31    format(i11)
  1954.     if(itd(it).eq.1)goto 33
  1955.     pa(12)='/'
  1956.     write(scr,31)itd(it)
  1957.     read(scr,36)(aa(kk),kk=1,11)
  1958. 36    format(22a1)
  1959.     i=12
  1960.     do 32 j=1,11
  1961.     if(aa(j).eq.' ')goto 32
  1962.     i=i+1
  1963.     if(i.ge.nst)call bomb(3,'print3  ')
  1964.     pa(i)=aa(j)
  1965. 32    continue
  1966. 33    j=mod(it+1-ia,1000)
  1967.     ir=23
  1968.     if(pa(ir).ne.' ')goto 41
  1969.     do 50 m=ir,4,-1
  1970.     if(pa(m).ne.' ')goto 51
  1971. 50    continue
  1972. 51    mm=min0(ir-m,8)
  1973.     do 55 i=ir,mm+1,-1
  1974.     pa(i)=pa(i-mm)
  1975. 55    pa(i-mm)=' '
  1976. 41    write(8,34)j,(pa(k),k=1,n)
  1977. 34    format(1x,i3,126a1,/,26x,106a1)
  1978. 40    continue
  1979.     return
  1980.     end
  1981.  
  1982.  
  1983.     subroutine prntap
  1984. $include: 'picom.for'
  1985.     do 1 n=nchpl,1,-1
  1986.     if(ap(n).ne.' ')goto 2
  1987. 1    continue
  1988.     n=1
  1989. 2    write(8,3)(ap(i),i=1,n)
  1990. 3    format(1x,127a1)
  1991.     return
  1992.     end
  1993.  
  1994.  
  1995.     subroutine read(ib)
  1996. $include: 'picom.for'
  1997. c  ib=0 on return if line has nonblanks
  1998. c  ib=1 on return if line is all blanks
  1999.     line=line+1
  2000.     ib=inpba(line)
  2001.     ie=inpea(line)
  2002.     il=ie-ib+1
  2003.     do 1 i=1,nchpl
  2004. 1    a(i)=' '
  2005.     call copyc(ainpl(ib),a,il)
  2006.     call copyc(a,ap,nchpl)
  2007.     i=1
  2008.     if(a(1).eq.'*')goto 8
  2009.     do 6 i=1,nchpl
  2010.     if(a(i).eq.';')goto 8
  2011. 6    continue
  2012.     goto 14
  2013. 8    do 10 j=i,nchpl
  2014. 10    a(j)=' '
  2015. 14    ib=0
  2016.     do 16 i=1,nchpl
  2017.     if(a(i).ne.' ')goto 18
  2018. 16    continue
  2019. 17    ib=1
  2020.     return
  2021. 18    if(i.eq.1)return
  2022.     do 20 j=i,nchpl
  2023.     a(j+1-i)=a(j)
  2024. 20    a(j)=' '
  2025.     return
  2026.     end
  2027.  
  2028.  
  2029.     subroutine readin
  2030. $include: 'picom.for'
  2031.     do 1 i=1,nlab
  2032. 1    linlab(i)=-1
  2033.     line=0
  2034.     lines=0
  2035.     labels=0
  2036.     iinpl=1
  2037. 5    read(7,6)a
  2038. 6    format(127a1)
  2039.     do 7 k=nchpl,1,-1
  2040.     if(a(k).ne.' ')goto 10
  2041. 7    continue
  2042.     k=1
  2043. 10    lines=lines+1
  2044.     lintot=lintot+1
  2045.     if(lines.gt.ninl)call bomb(4,'NINL    ')
  2046.     if(k.eq.1)goto 15
  2047.     call srchnb(1,k,m)
  2048.     if(a(m).eq.'p' .and. a(m+1).eq.'r' .and. a(m+2).eq.'o')goto 11
  2049.     if(a(m).ne.'l' .or. a(m+1).ne.'a' .or. a(m+2).ne.'b')goto 15
  2050. 11    call nxt(m+2,nchpl,' ',j1)
  2051.     call nxtch(j1,nchpl,j3)
  2052.     call nxt(j3,nchpl,' ',j4)
  2053.     call numl(j1,nchpl,j3,j4,j5)
  2054.     if(linlab(j5).ne.-1)write(8,12)(a(i),i=j3,j4)
  2055. 12    format(' duplicate label: ', 15a1)
  2056.     if(linlab(j5).ne.-1)call bomb(3,' readin1')
  2057.     linlab(labels)=lines
  2058. 15    inpba(lines)=iinpl
  2059.     inpea(lines)=iinpl+k-1
  2060.     if(inpea(lines).gt.ninch)call bomb(4,'NINCH   ')
  2061.     call copyc(a,ainpl(iinpl),k)
  2062.     iinpl=iinpl+k
  2063.     if(k.eq.1)goto 5
  2064.     if(a(m).eq.'f' .and. a(m+1).eq.'l' .and. a(m+2).eq.'u')goto 16
  2065.     if(a(m).ne.'e' .or. a(m+1).ne.'n' .or. a(m+2).ne.'d')goto 5
  2066.     if(a(m+3).ne.' ' .or. a(m+4).ne.'i' .or. a(m+5).ne.'n')goto 5
  2067. 16    write(8,20)lines,labels,iinpl
  2068. 20    format(' number of lines in input file=',i5,' labels='
  2069.      +      ,i2,'  characters=',i6)
  2070.     return
  2071.     end
  2072.  
  2073.  
  2074.     subroutine recip(j13,j4)
  2075. $include: 'picom.for'
  2076.     integer*4 in,id
  2077.     ia=isb(j13)
  2078.     iz=ise(j13)
  2079.     k=itb(ia)
  2080.     if(ivn(k).eq.2 .and. ivp(k).ne.0 .and. ia.ne.iz)then
  2081.         write(8,2)
  2082.         call bomb(3,'  recip1')
  2083. 2    format(' this simple program cannot process the denominator.',
  2084.      +     '  denominators should have 1 term without variable 1')
  2085.     endif
  2086.     ja=nt+1
  2087.     call copyt(ia,ia,1,1)
  2088.     in=itn(ja)
  2089.     id=itd(ja)
  2090.     itn(ja)=id
  2091.     if(in.lt.0)itn(ja)=-id
  2092.     itd(ja)=iabs(in)
  2093.     do 4 i=itb(ja),ite(ja)
  2094. 4    ivp(i)=-ivp(i)
  2095.     if(iz.eq.ia)call defsum(j4,ja,nt)
  2096.     if(iz.eq.ia)return
  2097.     if(maxe.gt.2)write(8,6)
  2098.     if(maxe.gt.2)call bomb(3,'  recip2')
  2099. 6    format(' this simple program cannot divide if maxe is',
  2100.      +      ' greater than 2')
  2101.     ib=ia+1
  2102.     k=itb(ib)
  2103.     if(ivn(k).ne.2 .or. ivp(k).lt.1)write(8,2)
  2104.     if(ivn(k).ne.2 .or. ivp(k).lt.1)call bomb(3,'  recip3')
  2105.     if(ivp(k).eq.2)goto 14
  2106.     ic=ib
  2107.     if(ib.eq.iz)goto 10
  2108. 8    k=itb(ic+1)
  2109.     if(ic.eq.iz .or. ivp(k).eq.2)goto 10
  2110.     ic=ic+1
  2111.     goto 8
  2112. 10    jb=nt+1
  2113.     call copyt(ib,ic,1,1)
  2114.     jc=nt
  2115.     if(ic.eq.iz)goto 12
  2116. 11    jd=nt+1
  2117.     call copyt(ic+1,iz,1,1)
  2118.     je=nt
  2119.     goto 16
  2120. 12    jd=nt+1
  2121.     je=jd
  2122.     itb(jd)=itetop+1
  2123.     call incnt
  2124.     call nite(itb(jd))
  2125.     k=itetop
  2126.     ivn(k)=2
  2127.     ivp(k)=2
  2128.     itn(jd)=0
  2129.     itd(jd)=1
  2130.     goto 16
  2131. 14    jb=nt+1
  2132.     jc=jb
  2133.     itb(jb)=itetop+1
  2134.     call incnt
  2135.     call nite(itb(jb))
  2136.     k=itetop
  2137.     ivn(k)=2
  2138.     ivp(k)=1
  2139.     itn(jb)=0
  2140.     itd(jb)=1
  2141.     ic=ib-1
  2142.     goto 11
  2143. 16    jf=nt+1
  2144.     call multn(ja,ja,jb,jc)
  2145.     jg=nt
  2146.     if(maxe.eq.1)goto 18
  2147.     jh=nt+1
  2148.     call multn(jf,jg,jf,jg)
  2149.     ji=nt
  2150.     jj=nt+1
  2151.     call multn(ja,ja,jd,je)
  2152.     jk=nt
  2153.     jl=nt+1
  2154.     call copyt(jf,jg,-1,1)
  2155.     call copyt(jh,ji,1,1)
  2156.     call copyt(jj,jk,-1,1)
  2157.     jm=nt
  2158.     call copyt(ja,ja,1,1)
  2159.     call multn(ja,ja,jl,jm)
  2160.     nt1=jm+1
  2161.     goto 24
  2162. 18    do 19 i=jf,jg
  2163. 19    itn(i)=-itn(i)
  2164.     call copyt(ja,ja,1,1)
  2165.     call multn(ja,ja,jf,jg)
  2166.     nt1=jg+1
  2167. 24    call defsum(j4,nt1,nt)
  2168.     call order(j4)
  2169.     return
  2170.     end
  2171.  
  2172.  
  2173.     subroutine swap(j4)
  2174. $include: 'picom.for'
  2175.     if=isb(j4)
  2176.     ig=ios(j4)
  2177.     if(if.ne.0)minns=min0(minns,ios(j4))
  2178.     ios(j4)=ios(1)
  2179.     iosi(ios(j4))=j4
  2180.     isb(j4)=isb(1)
  2181.     ise(j4)=ise(1)
  2182.     if(if.eq.0)goto 5
  2183.     ios(1)=ig
  2184.     iosi(ig)=1
  2185.     isb(1)=if
  2186.     ise(1)=if
  2187.     ite(if)=itb(if)
  2188.     call pack
  2189.     return
  2190. 5    call niltrm
  2191.     isb(1)=0
  2192.     call defsum(1,nt,nt)
  2193.     return
  2194.     end
  2195.  
  2196.  
  2197.     subroutine srch(j1,j2,h,j3)
  2198. $include: 'picom.for'
  2199.     character*1 h
  2200.     do 1 j3=j1,j2
  2201.     if(a(j3).eq.h)return
  2202. 1    continue
  2203.     j3=-1
  2204.     return
  2205.     end
  2206.  
  2207.  
  2208.     subroutine srchch(j1,j2,j3)
  2209. $include: 'picom.for'
  2210.     do 2 j3=j1,j2
  2211.     do 1 i=1,26
  2212.     if(a(j3).eq.ha(i))return
  2213.     if(a(j3).eq.hc(i))return
  2214. 1    continue
  2215. 2    continue
  2216.     j3=-1
  2217.     end
  2218.  
  2219.  
  2220.     subroutine srchnb(j1,j2,j3)
  2221. $include: 'picom.for'
  2222.     do 1 j3=j1,j2
  2223.     if(a(j3).ne.' ')return
  2224. 1    continue
  2225.     j3=-1
  2226.     return
  2227.     end
  2228.  
  2229.  
  2230.     subroutine sub
  2231. $include: 'picom.for'
  2232.     integer*4 inc4,idc4
  2233.     character*4 aj4,aj8
  2234.     call nxt(1,nchpl,' ',j1)
  2235.     call nxtnb(j1,nchpl,j2)
  2236.     call nxtch(j1,nchpl,j3)
  2237.     if(j2.ne.j3)goto 10
  2238.     call numv(j1,nchpl,j2,j3,j4)
  2239.     aj4='sum'
  2240.     if(isb(j4).eq.0)aj4='var'
  2241. 1    call nxt(j3+1,nchpl,'f',j5)
  2242.     if(a(j5+1).ne.'o')call bomb(3,' readin2')
  2243.     if(a(j5+2).ne.'r')call bomb(3,' readin4')
  2244.     if(a(j5+3).ne.' ')call bomb(3,' readin6')
  2245.     call numv(j5+3,nchpl,j6,j7,j8)
  2246.     aj8='var'
  2247.     if(isb(j8).gt.0)aj8='sum'
  2248.     call nxt(j7+1,nchpl,'i',j9)
  2249.     if(a(j9+1).ne.'n')call bomb(3,' readin8')
  2250.     if(a(j9+2).ne.' ')call bomb(3,'readin10')
  2251.     call numv(j9+2,nchpl,j10,j11,j12)
  2252.     minns=min0(minns,ios(j12))
  2253.     if(aj4.eq.'num')goto 35
  2254.     if(aj4.eq.'var')goto 30
  2255.     if(aj8.eq.'sum')goto 50
  2256.     call subs(j4,j8,j12)
  2257.     return
  2258. 10    call coef(j1,j4,inc4,idc4)
  2259.     inc=inc4
  2260.     idc=idc4
  2261.     aj4='num'
  2262.     j3=j4-1
  2263.     goto 1
  2264. 30    if(aj8.eq.'sum')goto 40
  2265.     call subv(j4,j8,j12)
  2266.     return
  2267. 35    if(aj8.eq.'sum')call bomb(3,'readin12')
  2268.     call subn(inc,idc,j8,j12)
  2269.     return
  2270. 40    call subr(j8,j4,j12)
  2271.     return
  2272. 50    a(71)='A'
  2273.     a(72)='?'
  2274.     a(73)='0'
  2275.     call numv(71,73,j15,j16,j17)
  2276.     call subr(j8,j17,j12)
  2277.     call subs(j4,j17,j12)
  2278.     return
  2279.     end
  2280.  
  2281.  
  2282.     subroutine subn(in,id,j8,j12)
  2283. $include: 'picom.for'
  2284.     integer*4 in4,id4
  2285.     if(isb(j12).eq.0)call bomb(3,'   subn1')
  2286.     if(isb(j8).ne.0)call bomb(3,'   subn2')
  2287.     do 8 it=isb(j12),ise(j12)
  2288.     do 5 i=itb(it),ite(it)
  2289.     if(ivn(i).ne.j8)goto 5
  2290.     in4=in
  2291.     id4=id
  2292.     in4=in4**ivp(i)
  2293.     id4=id4**ivp(i)
  2294.     call mulrat(itn(it),itd(it),in4,id4,itn(it),itd(it))
  2295.     ivp(i)=0
  2296.     if(ite(it).eq.itb(it))goto 5
  2297.     ite(it)=ite(it)-1
  2298.     if(i.gt.ite(it))goto 5
  2299.     do 1 j=i,ite(it)
  2300.     ivn(j)=ivn(j+1)
  2301. 1    ivp(j)=ivp(j+1)
  2302. 5    continue
  2303. 8    continue
  2304.     call order(j12)
  2305.     return
  2306.     end
  2307.  
  2308.  
  2309.     subroutine subr(ia,ib,ic)
  2310. c  replace pattern=first term in sum ia by variable ib in sum ic
  2311. $include: 'picom.for'
  2312.     integer*4 inj1,idj1
  2313.     if(isb(ia).eq.0)call bomb(3,'   subr1')
  2314.     if(isb(ib).ne.0)call bomb(3,'   subr2')
  2315.     nt1=nt+1
  2316.     j1=isb(ia)
  2317.     j2=itb(j1)
  2318.     j3=ite(j1)
  2319.     j4=isb(ic)
  2320.     j5=ise(ic)
  2321.     do 30 it=j4,j5
  2322.     j6=itb(it)
  2323.     j7=ite(it)
  2324.     minrat=30000
  2325.     do 5 j8=j2,j3
  2326.     jfp=ivn(j8)
  2327. 1    jfm=ivn(j6)
  2328.     if(jfm.gt.jfp)goto 22
  2329.     irat=ivp(j6)/ivp(j8)
  2330.     if(jfm.eq.jfp .and. irat.ge.1)goto 4
  2331.     jt=it+1-j4
  2332.     if(j2.eq.j3 .and. np.ge.2 .and. jfm.eq.jfp .and. irat.le.-1)
  2333.      +  write(8,2)jt,(pv(i,ic),i=1,nch),(pv(i,jfp),i=1,nch),
  2334.      +          (pv(i,ia),i=1,nch)
  2335. 2    format(' term',i4,' of ', 15a1,' has exponent of ', 15a1,
  2336.      +       ' with different sign from ', 15a1)
  2337.     j6=j6+1
  2338.     if(j6.le.j7)goto 1
  2339.     goto 22
  2340. 4    minrat=min0(minrat,ivp(j6)/ivp(j8))
  2341. 5    continue
  2342.     call incnt
  2343.     idj1=itd(j1)
  2344.     if(itn(j1).lt.0)idj1=-idj1
  2345.     inj1=iabs(itn(j1))
  2346.     call mulrat(itn(it),itd(it),idj1,inj1,itn(nt),itd(nt))
  2347.     itb(nt)=itetop+1
  2348.     j6=itb(it)
  2349.     k=itb(nt)
  2350.     do 20 j=j6,j7
  2351.     do 12 j8=j2,j3
  2352.     if(ivn(j8).eq.ivn(j))goto 16
  2353. 12    continue
  2354.     ivn(k)=ivn(j)
  2355.     ivp(k)=ivp(j)
  2356.     goto 20
  2357. 16    ivn(k)=ivn(j)
  2358.     ivp(k)=ivp(j)-minrat*ivp(j8)
  2359. 20    k=k+1
  2360.     ivn(k)=ib
  2361.     ivp(k)=minrat
  2362.     call nite(k)
  2363.     goto 30
  2364. 22    call copyt(it,it,1,1)
  2365. 30    continue
  2366.     call defsum(ic,nt1,nt)
  2367.     call order(ic)
  2368.     return
  2369.     end
  2370.  
  2371.  
  2372.     subroutine subs(in1,in2,in3)
  2373. c  substitute sum=in1 for variable=in2 in sum=in3
  2374. c  result goes in nt+1,...
  2375. $include: 'picom.for'
  2376.     integer*4 inp,idp
  2377.     common/nwt/iitb,iite,jitb,jite,ntb,if
  2378. c  dimension of ittbs ittes is max jp
  2379.     dimension ittbs(30),ittes(30)
  2380.     jp=1
  2381.     ia=in2
  2382.     ib=isb(in3)
  2383.     ic=ise(in3)
  2384.     id=isb(in1)
  2385.     ie=ise(in1)
  2386.     if(ib.eq.0 .or. id.eq.0)call bomb(3,'   subs1')
  2387.     ntb=nt+1
  2388.     lev1=0
  2389.     do 30 it=ib,ic
  2390.     iitb=itb(it)
  2391.     iite=ite(it)
  2392.     if(ivn(iitb).eq.2)lev1=ivp(iitb)
  2393.     if(lev1.gt.maxe)goto 31
  2394.     do 5 if=iitb,iite
  2395.     if(ivn(if).ne.ia)goto 5
  2396.     ifp=ivp(if)
  2397.     lev2=0
  2398.     if(ifp.eq.1)goto 10
  2399.     if(id.eq.ie)goto 70
  2400.     if(ifp.gt.0)goto 17
  2401.     itpr=it-ib+1
  2402.     if(np.ge.2)write(8,1)(pv(i,ia),i=1,nch),
  2403.      +       (pv(i,in3),i=1,nch),itpr
  2404. 1    format(2x, 15a1,' occurs to a negative power in ', 15a1
  2405.      +    , '  term ',i4)
  2406. 5    continue
  2407. c  here copy term it to nt+1
  2408.     if=0
  2409.     jitb=0
  2410.     call sub2(itn(it),itd(it))
  2411.     goto 30
  2412. c  here copy term it with substitution
  2413. 10    do 68 jt=id,ie
  2414.     jitb=itb(jt)
  2415.     jite=ite(jt)
  2416.     if(ivn(jitb).eq.2)lev2=ivp(jitb)
  2417.     if(lev1+lev2.gt.maxe)goto 30
  2418.     call mulrat(itn(it),itd(it),itn(jt),itd(jt),inp,idp)
  2419.     call sub2(inp,idp)
  2420. 68    continue
  2421.     goto 30
  2422. 70    if(ivn(itb(id)).eq.2)lev2=ivp(itb(id))*ifp
  2423.     if(lev1+lev2.gt.maxe)goto 30
  2424.     k=0
  2425.     do 71 i=itb(id),ite(id)
  2426.     k=k+1
  2427.     jvn(k)=ivn(i)
  2428. 71    jvp(k)=ivp(i)*ifp
  2429.     jitb=1
  2430.     jite=k
  2431.     if(ifp.lt.0)goto 72
  2432.     call mulrat(itn(it),itd(it),itn(id)**ifp,itd(id)**ifp,inp,idp)
  2433.     goto 73
  2434. 72    ii=-ifp
  2435.     call mulrat(itn(it),itd(it),itd(id)**ii,itn(id)**ii,inp,idp)
  2436. 73    call sub3(inp,idp)
  2437.     goto 30
  2438. c  here get in1**ifp
  2439. 17    jd=id
  2440.     je=ie
  2441.     ifpsav=ifp
  2442.     ntbsav=ntb
  2443.     iitbsv=iitb
  2444.     iitesv=iite
  2445.     ifsav=if
  2446.     ntbmu=nt
  2447.     itesav=itetop
  2448. c  test whether this power is in memory
  2449.     if(jp.eq.1)goto 80
  2450.     if(ifp.lt.jp)goto 27
  2451.     if(ittbs(jp).eq.0)goto 30
  2452.     if(ifp.eq.jp)goto 27
  2453.     jd=nt+1
  2454.     k=itetop
  2455.     do 86 i=ittbs(jp),ittes(jp)
  2456.     call incnt
  2457.     itn(nt)=jtn(i)
  2458.     itd(nt)=jtd(i)
  2459.     itb(nt)=k+1
  2460.     do 84 j=jtb(i),jte(i)
  2461.     k=k+1
  2462.     ivn(k)=jvn(j)
  2463. 84    ivp(k)=jvp(j)
  2464.     call nite(k)
  2465. 86    continue
  2466.     je=nt
  2467.     goto 18
  2468. 80    ks=0
  2469.     nts=0
  2470. 18    jpb=nt+1
  2471.     ntbtmu=nt
  2472.     call mult(id,ie,jd,je)
  2473.     if(nt.eq.ntbtmu)then
  2474.         jp=jp+1
  2475.         ittbs(jp)=0
  2476.         nt=ntbmu
  2477.         ntb=ntbsav
  2478.         itetop=itesav
  2479.         goto 30
  2480.         endif
  2481.     jpe=nt
  2482.     jp=jp+1
  2483.     if(jp.gt.30)call bomb(4,'  max jp')
  2484.     jd=jpb
  2485.     je=jpe
  2486.     ittbs(jp)=nts+1
  2487.     do 21 i=jpb,jpe
  2488.     nts=nts+1
  2489.     jtn(nts)=itn(i)
  2490.     jtd(nts)=itd(i)
  2491.     jtb(nts)=ks+1
  2492.     do 20 j=itb(i),ite(i)
  2493.     ks=ks+1
  2494.     jvn(ks)=ivn(j)
  2495. 20    jvp(ks)=ivp(j)
  2496.     jte(nts)=ks
  2497.     if(ks.gt.ntott)call bomb(4,'NTOTT   ')
  2498.     if(nts.gt.ntermt)call bomb(4,'NTERMT  ')
  2499. 21    continue
  2500.     ittes(jp)=nts
  2501.     if(jp.ne.ifp)goto 18
  2502.     nt=ntbmu
  2503.     itetop=itesav
  2504.     if=ifsav
  2505.     ifp=ifpsav
  2506.     iitb=iitbsv
  2507.     iite=iitesv
  2508.     ntb=ntbsav
  2509. 27    do 54 jt=ittbs(ifp),ittes(ifp)
  2510.     jitb=jtb(jt)
  2511.     jite=jte(jt)
  2512.     if(jvn(jitb).eq.2)lev2=jvp(jitb)
  2513.     if(lev1+lev2.gt.maxe)goto 30
  2514.     call mulrat(itn(it),itd(it),jtn(jt),jtd(jt),inp,idp)
  2515.     call sub3(inp,idp)
  2516. 54    continue
  2517. 30    continue
  2518. 31    call defsum(in3,ntb,nt)
  2519.     return
  2520.     end
  2521.  
  2522.  
  2523.     subroutine subv(j4,j8,j12)
  2524. $include: 'picom.for'
  2525.     if(isb(j12).eq.0)call bomb(3,'   subv1')
  2526.     if(isb(j4)+isb(j8).ne.0)call bomb(3,'   subv2')
  2527.     do 8 it=isb(j12),ise(j12)
  2528.     do 5 i=itb(it),ite(it)
  2529.     if(ivn(i).eq.j8)ivn(i)=j4
  2530. 5    continue
  2531. 8    continue
  2532.     call order(j12)
  2533.     return
  2534.     end
  2535.  
  2536.  
  2537.     subroutine sub2(in,id)
  2538. $include: 'picom.for'
  2539.     common/nwt/iitb,iite,jitb,jite,ntb,if
  2540.     integer*4 in,id
  2541.     if(in.eq.0)return
  2542.     call incnt
  2543.     k=itetop
  2544.     itb(nt)=k+1
  2545.     itn(nt)=in
  2546.     itd(nt)=id
  2547.     ii=iitb
  2548.     if(jitb.eq.0)goto 15
  2549.     jj=jitb
  2550. 11    if(ii.gt.iite)goto 13
  2551.     if(ii.eq.if)goto 12
  2552.     if(jj.gt.jite)goto 15
  2553.     if(ivn(ii).lt.ivn(jj))goto 16
  2554.     if(ivn(ii).eq.ivn(jj))goto 57
  2555.     if(ivp(jj).eq.0)goto 18
  2556.     k=k+1
  2557.     ivn(k)=ivn(jj)
  2558.     ivp(k)=ivp(jj)
  2559. 18    jj=jj+1
  2560.     goto 11
  2561. 12    ii=ii+1
  2562.     goto 11
  2563. 13    if(jj.gt.jite)goto 60
  2564.     if(ivp(jj).eq.0)goto 17
  2565.     k=k+1
  2566.     ivn(k)=ivn(jj)
  2567.     ivp(k)=ivp(jj)
  2568. 17    jj=jj+1
  2569.     goto 13
  2570. 14    if(ii.gt.iite)goto 60
  2571. 15    if(ivp(ii).eq.0)goto 9
  2572.     k=k+1
  2573.     ivn(k)=ivn(ii)
  2574.     ivp(k)=ivp(ii)
  2575. 9    ii=ii+1
  2576.     if(ii.eq.if)goto 9
  2577.     goto 14
  2578. 16    if(ivp(ii).eq.0)goto 12
  2579.     k=k+1
  2580.     ivn(k)=ivn(ii)
  2581.     ivp(k)=ivp(ii)
  2582.     ii=ii+1
  2583.     goto 11
  2584. 57    if(ivp(ii)+ivp(jj).eq.0)goto 58
  2585.     k=k+1
  2586.     ivn(k)=ivn(ii)
  2587.     ivp(k)=ivp(ii)+ivp(jj)
  2588. 58    ii=ii+1
  2589.     jj=jj+1
  2590.     goto 11
  2591. 60    if(k.gt.itetop)goto 56
  2592.     k=itetop+1
  2593.     ivn(k)=2
  2594.     ivp(k)=0
  2595. 56    itesav=itetop
  2596.     call nite(k)
  2597.     if(nt.eq.ntb)goto 68
  2598.     if(icomp(ntb,nt))31,33,34
  2599. 31    call tzap(nt,ntb)
  2600.     goto 68
  2601. 33    call addrat(itn(ntb),itd(ntb),itn(nt),itd(nt),itn(ntb),itd(ntb))
  2602.     nt=nt-1
  2603.     itetop=itesav
  2604.     goto 68
  2605. 34    if(icomp(nt-1,nt))62,61,68
  2606. 61    call addrat(itn(nt-1),itd(nt-1),itn(nt),itd(nt),
  2607.      ,     itn(nt-1),itd(nt-1))
  2608.     nt=nt-1
  2609.     itetop=itesav
  2610.     goto 68
  2611. 62    la=ntb
  2612.     lb=nt-1
  2613. 63    if(lb-la.eq.1)goto 64
  2614.     lc=(la+lb)/2
  2615.     if(icomp(lc,nt))65,66,67
  2616. 64    call tzap(nt,lb)
  2617.     goto 68
  2618. 65    lb=lc
  2619.     goto 63
  2620. 66    call addrat(itn(lc),itd(lc),itn(nt),itd(nt),itn(lc),itd(lc))
  2621.     nt=nt-1
  2622.     itetop=itesav
  2623.     goto 68
  2624. 67    la=lc
  2625.     goto 63
  2626. 68    return
  2627.     end
  2628.  
  2629.  
  2630.     subroutine sub3(in,id)
  2631. $include: 'picom.for'
  2632.     integer*4 in,id
  2633.     common/nwt/iitb,iite,jitb,jite,ntb,if
  2634.     if(in.eq.0)return
  2635.     call incnt
  2636.     k=itetop
  2637.     itb(nt)=k+1
  2638.     itn(nt)=in
  2639.     itd(nt)=id
  2640.     ii=iitb
  2641.     if(jitb.eq.0)goto 15
  2642.     jj=jitb
  2643. 11    if(ii.gt.iite)goto 13
  2644.     if(ii.eq.if)goto 12
  2645.     if(jj.gt.jite)goto 15
  2646.     if(ivn(ii).lt.jvn(jj))goto 16
  2647.     if(ivn(ii).eq.jvn(jj))goto 57
  2648.     if(jvp(jj).eq.0)goto 18
  2649.     k=k+1
  2650.     ivn(k)=jvn(jj)
  2651.     ivp(k)=jvp(jj)
  2652. 18    jj=jj+1
  2653.     goto 11
  2654. 12    ii=ii+1
  2655.     goto 11
  2656. 13    if(jj.gt.jite)goto 60
  2657.     if(jvp(jj).eq.0)goto 17
  2658.     k=k+1
  2659.     ivn(k)=jvn(jj)
  2660.     ivp(k)=jvp(jj)
  2661. 17    jj=jj+1
  2662.     goto 13
  2663. 14    if(ii.gt.iite)goto 60
  2664. 15    if(ivp(ii).eq.0)goto 9
  2665.     k=k+1
  2666.     ivn(k)=ivn(ii)
  2667.     ivp(k)=ivp(ii)
  2668. 9    ii=ii+1
  2669.     if(ii.eq.if)goto 9
  2670.     goto 14
  2671. 16    if(ivp(ii).eq.0)goto 12
  2672.     k=k+1
  2673.     ivn(k)=ivn(ii)
  2674.     ivp(k)=ivp(ii)
  2675.     ii=ii+1
  2676.     goto 11
  2677. 57    if(ivp(ii)+jvp(jj).eq.0)goto 58
  2678.     k=k+1
  2679.     ivn(k)=ivn(ii)
  2680.     ivp(k)=ivp(ii)+jvp(jj)
  2681. 58    ii=ii+1
  2682.     jj=jj+1
  2683.     goto 11
  2684. 60    if(k.gt.itetop)goto 56
  2685.     k=itetop+1
  2686.     ivn(k)=2
  2687.     ivp(k)=0
  2688. 56    itesav=itetop
  2689.     call nite(k)
  2690.     if(nt.eq.ntb)goto 68
  2691.     if(icomp(ntb,nt))31,33,34
  2692. 31    call tzap(nt,ntb)
  2693.     goto 68
  2694. 33    call addrat(itn(ntb),itd(ntb),itn(nt),itd(nt),itn(ntb),itd(ntb))
  2695.     nt=nt-1
  2696.     itetop=itesav
  2697.     goto 68
  2698. 34    if(icomp(nt-1,nt))62,61,68
  2699. 61    call addrat(itn(nt-1),itd(nt-1),itn(nt),itd(nt),
  2700.      ,     itn(nt-1),itd(nt-1))
  2701.     nt=nt-1
  2702.     itetop=itesav
  2703.     goto 68
  2704. 62    la=ntb
  2705.     lb=nt-1
  2706. 63    if(lb-la.eq.1)goto 64
  2707.     lc=(la+lb)/2
  2708.     if(icomp(lc,nt))65,66,67
  2709. 64    call tzap(nt,lb)
  2710.     goto 68
  2711. 65    lb=lc
  2712.     goto 63
  2713. 66    call addrat(itn(lc),itd(lc),itn(nt),itd(nt),itn(lc),itd(lc))
  2714.     nt=nt-1
  2715.     itetop=itesav
  2716.     goto 68
  2717. 67    la=lc
  2718.     goto 63
  2719. 68    return
  2720.     end
  2721.  
  2722.  
  2723.     subroutine switch(j,k,is)
  2724. $include: 'picom.for'
  2725.     do 2 it=isb(is),ise(is)
  2726.     do 1 i=itb(it),ite(it)
  2727.     if(ivn(i).eq.j)ivn(i)=k
  2728. 1    continue
  2729. 2    continue
  2730.     return
  2731.     end
  2732.  
  2733.  
  2734.  
  2735.     subroutine taylor
  2736. $include: 'picom.for'
  2737.     character*4 star
  2738.     integer*4 itf
  2739.     dimension ltay(mitay)
  2740.     data iorder,ifirst/1,0/
  2741.     if(idflg.ne.0)call dertab
  2742.     call nxt(2,nchpl,' ',j1)
  2743.     call srch(j1,nchpl,'=',j2)
  2744.     j4=j1
  2745.     if(j2.eq.-1)goto 5
  2746.     call getexp(j2+1,iorder,j3,0)
  2747.     if(iorder.lt.1)call bomb(3,' taylor1')
  2748.     return
  2749. 5    call nxt(2,nchpl,' ',j13)
  2750.     star='no'
  2751.     if(a(j13-1).eq.'*')star='yes'
  2752. 10    call nxt(j4,nchpl,'i',j3)
  2753.     if(j3.eq.-1)call bomb(3,' taylor2')
  2754.     if(a(j3-1).eq.' ' .and. a(j3+1).eq.'n' .and. a(j3+2).eq.' ')
  2755.      +       goto 12
  2756.     j4=j3+1
  2757.     if(j4.lt.nchpl)goto 10
  2758.     call bomb(3,' taylor3')
  2759. 12    do 14 j5=j3-1,j1,-1
  2760.     if(a(j5).ne.' ')goto 16
  2761. 14    continue
  2762.     call bomb(3,' taylor4')
  2763. 16    do 18 j6=j5,j1,-1
  2764.     if(a(j6).eq.' ')goto 20
  2765. 18    continue
  2766.     call bomb(3,' taylor5')
  2767. 20    call numv(j6,j5,j7,j8,j9)
  2768.     if(isb(j9).eq.0)call bomb(3,' taylor6')
  2769.     itesav=itetop
  2770.     call term(j3+2)
  2771.     itay=itetop-itesav
  2772.     if(itay.gt.mitay)call bomb(4,'MITAY   ')
  2773.     do 21 i=1,itay
  2774. 21    ltay(i)=ivn(itesav+i)
  2775.     nt=nt-1
  2776.     itetop=itesav
  2777.     call niltrm
  2778.     call defsum(1,nt,nt)
  2779.     if(ifirst.ne.0)goto 22
  2780.     a(1)='A'
  2781.     a(2)='?'
  2782.     a(3)='0'
  2783.     call numv(1,3,j14,j15,j10)
  2784.     a(3)='1'
  2785.     call numv(1,3,j14,j15,j11)
  2786.     a(3)='2'
  2787.     call numv(1,3,j14,j15,j12)
  2788. 22    ifirst=1
  2789.     nv1=nv+1
  2790.     do 34 it=isb(j9),ise(j9)
  2791.     call copyt(it,it,1,1)
  2792.     call defsum(j10,nt,nt)
  2793.     do 32 ivv=1,itay
  2794.     iv=ltay(ivv)
  2795.     if(star.eq.'yes')call switch(iv,nv1,j10)
  2796.     call copys(j10,j11)
  2797.     call subn(0,1,iv,j10)
  2798.     itf=1
  2799.     do 30 io=1,iorder
  2800.     call deriv(j11,iv,j12)
  2801.     call outqu(j12)
  2802.     if(itn(isb(j12)).eq.0)goto 31
  2803.     call copys(j12,j11)
  2804.     call subn(0,1,iv,j12)
  2805.     itf=itf*io
  2806.     call incnt
  2807.     k=itetop+1
  2808.     itn(nt)=1
  2809.     itd(nt)=itf
  2810.     itb(nt)=k
  2811.     ite(nt)=k
  2812.     ivn(k)=iv
  2813.     ivp(k)=io
  2814.     call nite(k)
  2815.     nt1=nt+1
  2816.     call mult(nt,nt,isb(j12),ise(j12))
  2817.     call defsum(j12,nt1,nt)
  2818.     call pack
  2819. 30    call adds(j12,j10,1,1,j10)
  2820. 31    if(star.eq.'yes')call switch(nv1,iv,j10)
  2821. 32    continue
  2822. 34    call adds(j10,1,1,1,1)
  2823.     call swap(j9)
  2824.     call order(j9)
  2825.     call delsum(j10)
  2826.     return
  2827.     end
  2828.  
  2829.  
  2830.     subroutine term(ist)
  2831. c  each term is on a separate line, free format with coefficient first
  2832. $include: 'picom.for'
  2833.     integer*4 inc,idc
  2834.     call incnt
  2835.     call coef(ist,j1,inc,idc)
  2836.     itn(nt)=inc
  2837.     itd(nt)=idc
  2838.     itb(nt)=itetop+1
  2839.     if(j1.eq.-1)then
  2840.         k=itb(nt)
  2841.         ivn(k)=2
  2842.         ivp(k)=0
  2843.         goto 10
  2844.         endif
  2845.     k=itetop
  2846. 5    isgn=1
  2847.     if(a(j1).eq.'*')j1=j1+1
  2848.     if(a(j1).eq.'/')isgn=-1
  2849.     if(a(j1).eq.'/')j1=j1+1
  2850.     call numv(j1,nchpl,j2,j3,j4)
  2851.     ie=1
  2852.     j5=-1
  2853.     if(j3.eq.nchpl-1 .and. a(nchpl).ne.' ')call bomb(3,'   term1')
  2854.     if(j3.le.nchpl-2)call getexp(j3+1,ie,j5,1)
  2855.     k=k+1
  2856.     ivn(k)=j4
  2857.     ivp(k)=ie*isgn
  2858.     j1=j5
  2859.     if(j5.ne.-1)goto 5
  2860. 10    call nite(k)
  2861.     return
  2862.     end
  2863.  
  2864.  
  2865.     subroutine thderv
  2866. $include: 'picom.for'
  2867.     idflg=1
  2868.     call nxt(2,nchpl,'o',j1)
  2869.     call nxt(j1,nchpl,' ',j2)
  2870.     call numv(j2,nchpl,j3,j4,j5)
  2871.     call nxt(j4+1,nchpl,'w',j6)
  2872.     call nxt(j6,nchpl,' ',j7)
  2873.     call numv(j7,nchpl,j8,j9,j10)
  2874.     call srch(j9+1,nchpl,'=',j11)
  2875.     if(j11.ne.-1)goto 1
  2876.     call srch(j9+1,nchpl,'i',j12)
  2877.     if(j12.eq.-1)call bomb(3,'thderv1')
  2878.     call nxt(j12,nchpl,' ',j11)
  2879. 1    call srchnb(j11+1,nchpl,j12)
  2880.     if(j12.eq.-1)call bomb(3,'thderv2')
  2881.     ig=0
  2882.     if(a(j12).eq.'0')ig=1
  2883.     if(ig.eq.1)goto 8
  2884.     call numv(j11+1,nchpl,j12,j13,j14)
  2885. 8    if=0
  2886.     do 10 i=1,nd
  2887.     if(ider(1,i).eq.j5 .and. ider(2,i).eq.j10)if=1
  2888.     if(ider(1,i).eq.j5 .and. ider(2,i).eq.j10)in=i
  2889. 10    continue
  2890.     if(if.eq.0 .and. ig.eq.1)goto 50
  2891.     if(if.eq.0)goto 14
  2892.     ja=ider(3,in)
  2893.     jb=ipv(ja)
  2894.     if(np.ge.2)write(8,11)(pv(i,ja),i=1,jb)
  2895. 11    format(' REDEFINITION: previously the  derivative was = ',15a1)
  2896.     if(ig.eq.0)goto 18
  2897.     nd=nd-1
  2898.     if(in.gt.nd)goto 50
  2899.     do 28 ih=in,nd
  2900.     ider(1,ih)=ider(1,ih+1)
  2901.     ider(2,ih)=ider(2,ih+1)
  2902. 28    ider(3,ih)=ider(3,ih+1)
  2903.     goto 50
  2904. 14    nd=nd+1
  2905.     in=nd
  2906.     if(nd.gt.nder)call bomb(4,'NDER    ')
  2907.     ider(1,in)=j5
  2908.     ider(2,in)=j10
  2909. 18    ider(3,in)=j14
  2910. 50    return
  2911.     end
  2912.  
  2913.  
  2914.     subroutine tzap(ka,kb)
  2915. c  zap term at ka to position kb
  2916. $include: 'picom.for'
  2917.     integer*4 itnka,itdka
  2918.     if(ka.eq.kb)return
  2919.     itnka=itn(ka)
  2920.     itdka=itd(ka)
  2921.     itbka=itb(ka)
  2922.     iteka=ite(ka)
  2923.     if(ka.gt.kb)goto 5
  2924.     do 1 i=ka,kb-1
  2925.     itn(i)=itn(i+1)
  2926.     itd(i)=itd(i+1)
  2927.     itb(i)=itb(i+1)
  2928. 1    ite(i)=ite(i+1)
  2929.     goto 8
  2930. 5    do 6 i=ka,kb+1,-1
  2931.     itn(i)=itn(i-1)
  2932.     itd(i)=itd(i-1)
  2933.     itb(i)=itb(i-1)
  2934. 6    ite(i)=ite(i-1)
  2935. 8    itn(kb)=itnka
  2936.     itd(kb)=itdka
  2937.     itb(kb)=itbka
  2938.     ite(kb)=iteka
  2939.     return
  2940.     end
  2941.  
  2942.  
  2943.  
  2944.     subroutine zsub(j3)
  2945. $include: 'picom.for'
  2946. c this is available
  2947.     return
  2948.     end
  2949.